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Abstract.  The  atmospheric  flow  generated  in  response  to  horizon¬ 
tal  gradients  in  the  surface  temperature  or  heatflux  is  investi¬ 
gated.  A  linear  model  is  described  that  can  treat  the  problem  for 
the  case  of  small  surface  temperature  gradients,  including  the 
approximate  effect  of  a  large-scale  flow  and  the  effect  of  a  ver¬ 
tical  structure  in  the  heat  and  momentum  diffusivity,  and  in  the 
stratification.  The  second  part  of  this  report  describes  a  nu¬ 
merical  mesoscale  model,  and  results  from  the  application  of 
this  model  to  the  simulation  of  the  sea-breeze. ^ 
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APPENDIX 
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The  sun  also  ariseth,  and  the  sun 
goeth  down,  and  hasteth  to  his 
place  where  he  arose. 

The  wind  goeth  toward  the  south 
and  turneth  about  unto  the  northy; 
it  whirleth  about  continually,  and 
the  wind  returneth  again  according 
to  his  circuits. 


Ecclesiastes  1:5,6 


1.  INTRODUCTION 


The  present  report  deals  with  some  aspects  of  the  atmospheric 
response  to  time-varying  horizontal  gradients  in  the  surface 
temperature  and  heatflux  on  the  mesoscale.  The  interest  is  con¬ 
centrated  on  the  modelling  of  the  sea-breeze  circulation  in  its 
"pure"  form  as  the  two-dimensional  flow  across  a  coastline  in 
response  to  the  periodic  heating  of  the  land. 

In  the  first  part  an  analytical  technique  is  employed,  which 
enables  the  study  of  different  scales  of  the  mean  flow.  The 
method  is  based  on  the  linearization  of  the  dynamical  equations, 
and  the  neglect  of  any  feedback  from  the  mean  flow  to  the  strati¬ 
fication  and  turbulence  structure  (constant  diffusivity  and 
stratification).  The  model  is  generalized  to  include  the  effect 
of  vertical  variation  of  diffusivity  and  stratification. 

The  second  part  describes  a  numerical  model  in  which  a  more 
sophisticated  treatment  of  the  turbulent  transport  is  employed. 

A  large  number  of  theoretical  studies  of  the  sea-breeze  phenom¬ 
enon  have  appeared  in  the  literature.  Most  of  these  are  based  on 
numerical  modelling  (Estoque  1961, Neumann  and  Mahrer  1971,Pielke 


1974,  and  many  others)  but  most  of  the  results  presented  in  these 
studies  relate  to  only  gross  features  of  the  flow  which  in  many 
cases  can  be  inferred  from  an  analysis  of  the  governing  equations 
using  analytical  techniques  (Walsh, 1974 ) .  The  accurate  modelling 
of  the  development  of  small-scale  structure  in  the  sea-breeze 
flow,  and  the  interaction  with  the  turbulence  field  is  however 
possible  only  by  using  a  numerical  technique.  In  recent  years 
the  interest  in  the  detailed  structure  of  the  flow  and  the  tur¬ 
bulence  fields  has  been  promoted  by  concern  for  the  implications 
to  dispersion  of  pollutants  in  the  coastal  zone  (Keen  and  Lyons 
1978,  Anthes  1978),  and  the  level  of  sophistication  of  oper¬ 
ational  dispersion  models,  in  particular  "puff"  models  (Mikkelsen 
et.  al.  1981),  has  allowed  for  a  direct  use  in  such  models  of 
the  detailed  windfields  and  turbulence  parameters  from  the  flow 
models.  The  model  presented  here  has  potential  for  such  practical 
use,  but  at  present  it  has  been  applied  only  for  simulations 
using  synthetic  "typical"  input  data. 

The  physical  situation  modelled  in  both  chapters,  except  for 
Section  2.9,  is  a  straight  coastline  dividing  a  heated  semiin¬ 
finite  land  area  and  a  sea  where  the  temperature  is  assumed 
constant.  The  coordinate  system  is  defined  as  a  right-handed 
one  with  the  x-axis  pointing  inland,  the  y-axis  along  the  coast¬ 
line,  and  z-axis  vertical.  The  flow  is  assumed  to  be  describable 
as  two-dimensional  with  no  variation  along  the  y-axis.  The 
components  of  the  wind  vector  u,  v,  and  w,  are,  as  usual,  de¬ 
fined  as  the  components  along  the  x-,  y-,  and  z-axis,  respect¬ 
ively. 
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2.  LINEARIZED  PERIODICALLY  PORCED  MODEL 


2.1.  The  scaled  basic  equations 


The  equations  of  motion  assuming  the  Boussinesq  approximation 
can  be  written  as: 


3$  1  3 

+  $  •  +  fit  x  J  «  _ Vp  +  t 

TE  ps  Tz 


36'  3$ 

„  „  +  $  •  V 9 1  +  wr  *  __  +  Q 

TfT  Tz 

(2.1) 

v3  •  v  -  0 


3p  0' 

Tz  "  9Pa  r  ’ 


Where  the  hydrostatic  approximation  has  been  also  assumed  and 
the  horizontal  diffusion  of  momentum  and  heat  have  been  neglected. 
The  hydrostatic  approximation  is  justified  if  the  aspect  ratio 
of  the  flow  can  be  assumed  to  be  much  less  than  one:  We  are  here 
considering  flow  driven  by  surface  heating  and  therefore  the 
pressure  gradient  is  the  important  driving  force  for  the  flow. 

We  can  therefore  assume  the  scaling: 


p  pU 

r  w 


(2.2) 


from  the  equation  for  the  horizontal  velocity,  where  L  is  a 
characteristic  length,  T  characteristic  time,  and  0  velocity. 
Prom  the  equations  of  continuity  and  vertical  velocity  we  can 
write: 
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W  HO  H2  p 

T  L  T  l?  H 


(2.3) 


Thus  the  magnitude  of  the  time  change  of  w  is  of  the  order 
(H/L)2  p-1  times  the  magnitude  of  the  vertical  pressure  gradi¬ 
ent.  Also  the  omission  of  the  horizontal  diffusion  of  momentum 
can  be  justified  assuming  H  <<  L.  In  a  viscous  fluid  we  have 
that  the  horizontal  diffusion  of  momentum  is  of  the  order  vu/L2, 
where  v  is  the  kinematic  viscosity  of  the  fluid;  the  vertical 
diffusion  term  is  similarly  of  the  order  vU/H2,  and  thus  the 
ratio  of  the  two  terms  is  (H/L)2.  In  the  case  of  the  diffusion 
by  turbulence  in  the  atmospheric  boundary  layer,  the  fluxes  of 
momentum  and  heat  are  typically  of  the  same  order  of  magnitude 
in  both  the  horizontal  and  vertical,  and  therefore  the  ratio  of 
the  two  terms  in  the  flux-divergence  terms  in  the  equation  of 
motion  and  in  the  thermodynamic  equation  is  of  the  order  (H/L). 


For  small  values  of  the  amplitude  of  the  thermal  forcing  the 
amplitude  of  the  flow  perturbation  will  be  small  and  the  non¬ 
linear  advection  terms  can  be  neglected.  The  equations  of  motion 
in  the  case  of  vanishing  basic  state  flow  becomes  in  this  case 
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where  the  further  simplification  has  been  made  that  the  vertical 
diffusion  terms  can  be  parameterized  employing  a  constant  dif- 
f usivity. 

We  assume  the  surface  forcing  to  be  harmonic  in  time  with  fre¬ 
quency  o>.  Equations  (2.4)  can  then  be  nondimensionalized  by  the 
introduction  of  the  following  scaling: 

■0  =  0*8 
p  =  Hg  ps  p 

g 

u  =  _  u 
N 


where  H  is  the  Brunt-V&is81M  frequency  corresponding  to  the  basic 
state  stratification,  and  variables  with  a  denote  dimensional 
quantities.  With  this  scaling  the  equations  of  motion  can  be 
written  as 
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where 

fg  ■  f/w 


I 
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(2.6) 
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2.2.  Method  of  solution 


In  the  following  we  specifically  consider  the  case  of  flow  driven 
by  diurnal  variation  of  surface  temperature  difference  across  a 
straight  coastline,  in  which  case  w  *  fl  *  angular  frequency  of 
the  earths  rotation;  thus  with  f  -  2  n  sin  ♦  we  have  f8  * 

2  sin  ♦  •  We  wish  to  obtain  periodic  solutions  to  this  system  of 
equations  with  the  following  boundary  conditions: 

at  z  ■  o  »0  :  u  ■  v  *  w  *  0 

0  for  5  <  0 

e  -  (2.7) 

cost  for  5  >  0 

at  t,  n  *  •  :  u,  v,  w,  p,  9  ♦  0 
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It  is  convenient  to  write  the  fields  in  terms  of  Fourier  series 
as  follows 


1  "  cos  k  5 

u(C,n,x)  *  _  /  u(k,n)  _  dk  •  exp(-ix)  (2.8) 

"  o  k 

and  similarly  for  v,  also 

1  "  sin  k? 

w(?,n,x)  ~  —  !  w(k,n)  — , -  dk  •  exp(-ix)  (2.9) 

»  o  K 

and  for  p  as  well.  The  boundary  condition  at  n  =  0  for  0  means 
that  6  must  be  expanded  as 


/ 1  1  “ 

9(5, nfx)  =  _  0o(n)  +  -  / 

\  2  *  o 


sin  k5 


B(k, n)dkjexp( -i x) 


(2.10) 


with  0  (k,0)  =  0O( 0 )  =  1.  Because  of  the  linearity  of  the  system 
it  can  be  solved  for  each  Fourier  mode  independently.  The  first 
term  on  the  right  in  Eq.  (2.10)  corresponds  to  the  heating  of 
the  surface  harmonically  in  time  and  independently  of  t.  This 
forcing  is  easily  seen  to  induce  no  flow,  but  to  give  contri¬ 
butions  to  the  0  and  p  fields  only.  When  only  this  forcing  term 
is  considered  the  0  equation  simplifies  to 


(‘ 


32  \ 
1  +  7^2/ 


e0  =  o 


(2.11) 


then 


e0(n)  =  exp^-  ^  n)cos(j-  n  -  x^ 


v/2 


(2.12) 


This  is  the  classical  Stokes  solution  to  the  problem  of  a  viscous 
fluid  subject  to  harmonic  forcing  at  the  surface.  The  equation 
governing  the  functions  u(k,  n)»  v(k,  n)  etc.  is  obtained  from 
the  set  of  Eqs.  (2.6): 
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For  small  values  of  fs,  corresponding  to  #  »  f  rotation  can  be 
neglected  and  Eq.  (2.14)  reduces  to 


36  34  a2 

7^*  2i  7^' 


(2.17) 


and  Eq.  (2.16)  becomes  in  this  case 


(i+a2)  a  ±  k  =  0 


(2.18) 


Here  we  are  interested  only  in  the  case  where  u>  =  1/(24  h),  and 
therefore  rotation  can  only  be  neglected  if  we  have  4>  £  10  deg. 

It  should  be  noted,  however,  that  if  the  surface  boundary  con¬ 
ditions  are  changed  so  that  the  horizontal  length  scale  no  longer 
is  determined  by  L  in  Eq.  (2.5),  but  rather  induced  externally  as 
l  with  l  <<  L,  then  only  large  values  of  k  in  Eq.  (2.14)  need 
be  considered.  In  this  case  Eq.  (2.14)  can  be  satisfied  only  if 
the  higher-order  terms  and  k2  are  of  comparable  magnitude  and 
thus  the  rotation  term  can  be  neglected.  It  is  obvious  that  the 
flow  structure  is  determined  in  this  case  solely  by  the  par¬ 
ameter  t/L.  A  discussion  of  the  flow  structure  in  this  case  is 
given  in  Kimura  and  Eguchi  (1978),  where  the  governing  parameter 
is  defined  as  n  =  (t/L)2/3.  They  neglect  rotation  even  in  the 
case  n  =  1,  however.  In  Eq.  (2.5)  we  have  introduced  the  time 
scale  as  w  ;  alternatively,  we  could  have  used  f-^.  The  equation 
for  the  perturbation  pressure  in  this  case  where  we  have  an  ex¬ 
ternally  specified  horizontal  length  scale  l  can  be  obtained 
from  Eq.  (2.14)  by  the  substitution  n'=  /w/T  n  yielding: 


36 

/LTTTi 


+  2iw 


+  (  1  —  li>4  ) 


i2 


p  =  0  (2.19) 


(2.20) 


In  this  case  the  vertical  scale  is  the  Ekman  depth  H  =  /K/f , 
n'  »  Z/H.  Eq.  (19)  makes  it  possible  to  identify  a  number  of 
prototype  boundary  layers  by  requiring  that  different  combi- 
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nations  of  two  terms  be  dominant.  This  analysis  is  performed  in 
Park  and  Mahrt  (1979),  where  it  is  argued  that  the  stratifi¬ 
cation  paramater  m  is  much  more  sensitive  to  the  choice  of  hori¬ 
zontal  length  scale  t  than  to  the  stratif ication.  This  latter 
term  is  given  by  N  because  the  diffusivity  K  should  not  be  chosen 
independently  of  N,  since  turbulence  is  suppressed  when  N  is 
large.  The  choice  of  i  >>  L  or  equivalently  m  <<  should,  how¬ 
ever,  be  expected  to  lead  to  incorrect  scaling  of  the  flow  par¬ 
ameters  because  even  though  a  typical  length  scale  characterizing 
the  heating/cooling  area  is  t,  then  the  "internal"  flow  scale  L 
is  the  proper  length  scale.  The  conclusion  is  thus  that  m  in 
actual  flows  is  limited  both  downwards  and  upwards. 

In  the  following  we  return  to  our  original  formulation,  and  are 
therefore  interested  in  the  solution  of  the  system  of  equations 
(2.16). 

The  solution  of  the  cubic  equation  is  easily  found  by  using 
Newton-Raphson  iteration  to  obtain  one  root  B^,  and  the  two  re¬ 
maining  roots  from  the  quadratic  obtained  when  (  B—  B -| )  is  factored 
out.  the  solution  for  the  a^'s  depends  on  the  wavenumber  k  and 
the  scaled  Coriolis  parameter  fs.  Figure  2.1  shows  the  real  and 
imaginary  part  of  the  three  values  of  a  from  the  solution  of  the 
cubic  as  function  of  k.  The  value  of  fs  is  1.5  corresponding  to 
diurnal  forcing  at  middle  latitudes  (♦  =  50°).  The  as/mptotic 
solution  corresponding  to  the  limit  k  +  »  is  given  as 

=  k’/3  exp(i(j  +  1)  I)  ,  j  *  1,2,3  (2.21) 

This  expression  is  a  fair  approximation  even  at  moderate  values 
of  k.  The  fourth  root  in  Eq.  (2.16)  is  independent  of  k  and 
given  as  a4  =  -  /2/2 { 1  —  i ) .  As  for  the  pressure  the  coefficients 
in  the  Fourier  integrals  Eqs.  (2.9)  and  (2.10)  can  be  written  as 
a  sum  of  four  exponentials  as  Eq.  (2.15).  The  coefficients  are 
related  through  the  following  relations  obtained  from  substi¬ 
tution  of  Eq.  (2.15)  into  Eqs.  (2.13): 


Fig.  2.1.  The  absolute  value  of  the  real  part  ( _ ) ,  and  the 

imaginary  part  ( - )  of  the  roots  a  of  the  characteristic 

equation  (2.16). 
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The  solution  for  pj  is  obtained  by  application  of  the  boundary 
conditions  Eq.  (2.7). 


16 


2.3.  Limit  of  large  wavenumbers 
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In  the  asymptotic  limit  of  k  ♦  •  we  can  use  aj  from  Eg.  (2.21) 
and  the  system  of  Eqs.  (2.22)  can  be  approximated  by  the  following 
system : 

l  exp^i(j  +  1)  ^pj  =  k~1/3 

[  exp^-  i(j  +  1)  =0  j  =  1,2,3  (2.23) 

l  exp  ( i ( j  +  1)*)pj  =  0 

and 


P4  x  “  f|(Pl  +  P2  +  • 


The  solution  is 


P4 


k-1/3 
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and 


u4  =  0 


The  cross  isobar  flow  field  u(k,  n)  is  therefore  given  as 
(Fig.  2.2)  for  large  k: 


u(k,n)  =  l  uj  exp(ajn) 


1  /T  /  1  \  cTT  "\ 

=  -  exp(-X)  +  exp^-  ,  Xj  cos'-  X 


(2.24) 


with 


x  =  kV3  n 


The  appropriate  dimensional  scale  height  is  thus  dependent  on 
the  horizontal  length  scale.  Denoting  this  by  *  we  have  k  =  L/t 
and  the  scale  height  therefore  becomes 


h 


(r)V3 


(2.25) 


In  Park  and  Mahrt  (1979)  the  dynamics  in  this  limit  of  small  * 
is  termed  E1/3-layer  dynamics.  In  Mahrt  (1979)  it  is  noted  that 
this  dynamics  may  describe  actual  small-scale  nocturnal  flows 
even  though  the  initial  conditions  are  complicated  since  the 
flow  structure  is  determined  only  by  the  balance  of  the  pressure 
gradient  and  momentum  diffusion.  This  is  brought  out  in  Eq. 
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Fig.  2.2.  The  Cross  isobar  component  u(k,n)  in  the  limit  of 
large  wavenumbers  k  as  function  of  k'/^n  f  where  n  is  the 
scaled  vertical  coordinate. 


(2.24),  showing  u(k,  n)  to  be  real,  and  therefore  the  cross 
isobar  flow  to  be  in  phase  with  the  temperature  perturbation. 


2.4.  Solution  for  small  wavenumbers 

In  the  limit  of  k2  <<  1  a  simple  solution  can  again  be  obtained; 
We  have  for  the  cubic  Eq.  (2.16)  the  approximate  roots 


fs  “  1 


+  0(k4) 


i(1  *  f„)  ± 


2f„  (1?fs) 


+  0(k4) 


(2.26) 


64  *  -  i  . 
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To  simplify  we  will  assume  fs  >  1  in  the  following:  Then 


a,  =  -  a  b  k  +  0(kJ) 


1  sj  , 

—  ( 1  +  i )  +  0  ( k^ ) 
b  2 


1  ST  , 

1  -  i)  +  0(k2) 

3  £ 


(2.27) 


a*  = 


Si  , 

(1  -  x, 


with 


a  =  ( f s  +  1)"1/2  and  b  =  ( f s  -  1)‘1/2 


Eqs.  (2.22)  and  the  surface  boundary  conditions.  Eqs.  (2.7) 
lead  to  the  following  linear  equations  for  the  expansion  coef¬ 
ficient  for  pj  (Eq.  2.15)  when  the  above  relations  are  inserted 
and  terms  only  up  to  first  order  in  wavenumber  k  are  retained: 
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Using  Eq.  (2.22)  the  leading  terms  in  the  cross  isobar  flow 
field  for  k  <<  1  becomes 
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retaining  only  the  first-order  terms  we  see  that  the  contri¬ 
bution  in  the  Fourier  integral  Eq.  (2.8)  is: 
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(2.31 ) 


The  U|-term  is  of  order  k^  and  can  be  neglected  for  small 
heights  n,  however,  since  the  root  a1  is  proportional  to  k 
(Eq.  2.27)  the  depth  of  the  disturbance  from  this  component 
increases  as  k  tends  to  zero.  For  decreasing  values  of  k,  there¬ 
fore,  part  of  the  return  flow  will  consist  of  an  increasingly 
deep  layer  of  very  weak  flow.  The  contribution  to  the  total 
profile  from  the  first  term  takes  the  form: 
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where 
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The  expression  (2.31)  can  be  rewritten  as 


,  C08  iv  £  , 

Re{u(k,n)  . .  ....  exp(-  It)  }  »  cos(k  £)  •  A(  n)  cos(  t+*(  n) ) 

k 

(2.33) 

In  Fig.  2.3  the  amplitude  A( n)  is  plotted  for  different  values 
of  the  scaled  Coriolis  parameter  f8.  The  phase  function  ♦( n)  is 
shown  in  Fig.  2.4.  For  large  values  of  fs  the  flow  structure 
becomes  asymptotically  independent  of  w,  as  should  be  expected, 
corresponding  to  a  quasi-steady  state.  The  proper  depth  scale 
becomes  the  Bkman  depth,  as  can  be  seen  from  Eq.  (2.31)  with 
f8  »  1 
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where 

"-IF— (rf/2 

Thus,  for  fs  >>  1  we  have  $(  n)  =  -  if/ 4  up  to  the  height  n  = 
n  •  (1/2fs)-1/2  where  it  becomes  equal  to  3n/4.  This  stepwise 
behaviour  is  seen  to  be  a  fair  approximation  to  #  even  at 
fs  =  2  (Fig.  2.4).  The  inflow  layer  becomes  deeper  as  fs  de¬ 
creases  towards  1 . 


Fig.  2.3.  The  amplitude  A( n)  for  the  u-component  at  small 
wavenumbers  (Eq.  2.33)  for  different  values  of  the  scaled 
Coriolis  parameter. 
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Fig.  2.4.  The  Phase  angle  ♦( n)  for  the  u-component  at  small 
wavenumbers  (Eq.  2.33)  for  different  values  of  the  scaled 
Coriolis  parameter. 


2.5.  Land-sea  breeze  circulation 


By  inspection  the  dynamics  of  the  small  scales  are  seen  to  be 
characterized  by  a  balance  between  momentum  diffusion  and  pres¬ 
sure  gradient  force  in  the  u-momentum  equation,  and  between  tem¬ 
perature  diffusion  and  adiabatic  cooling  in  the  thermodynamic 
equation.  By  contrast  the  large  scales  are  found  to  be  charac¬ 
terized  by  a  balance  between  the  harmonic  term,  Coriolis  force, 
and  the  vertical  momentum  diffusion  for  the  u  and  v  components, 
whereas  in  the  thermodynamic  equation  the  adiabatic  cooling  is 
unimportant  and  we  have  the  Stokes  balance  between  forcing  and 
vertical  diffusion.  Well  above  the  surface,  however,  the  dif¬ 
fusion  terms  are  unimportant  and  the  principal  momentum  balance 
is  between  the  harmonic  terms,  Coriolis  term,  and  the  pressure 
gradient,  whereas  the  thermodynamic  equation  here  simplifies  to 
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a  balance  between  the  harmonic  term  and  the  adiabatic  cooling 
term. 

For  intermediate  wavenumbers  all  terms  are  important  and  even 
though  the  cubic  characteristic  equation  can  be  solved  to  yield 
explicit  roots  for  any  wavenumber,  the  final  solutions  obtained 
after  the  solution  of  the  linear  equations  for  the  boundary 
conditions  are  so  complicated  that  explicit  expressions  are  of 
little  use.  The  amplitude  and  phase  for  different  wavenumbers 
are  shown  in  Fig.  2.5.  Also  shown  is  the  amplitude  and  phase 
functions  for  the  integrated  profile  at  the  coast  obtained  by 
numerically  integrating  Eq.  (2.8)  using  Simpsons  rule.  The 
upper  limit  was  set  to  k  =  100  and  stepsize  in  the  integration 
was  Ak  =  0.1.  The  profiles  corresponding  to  individual  wave- 
numbers  are  characterized  by  a  nearly  constant  "inflow"  layer 
where  the  phase  is  essentially  constant  with  height?  at  the 
top  of  this  layer  the  amplitude  is  close  to  zero  and  the  phase 
changes  180  degrees  over  a  rather  thin  layer.  The  profiles 


Q  . .  1  *— 

0  0.1  0.2 


up))  AMPLITUDE 


Fig.  2.5a.  The  amplitude  of  the  modes  of  the  u-component 
for  different  wavenumbers.  The  dotted  curve  gives  the 
amplitude  of  the  integrated  u-profile. 
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plotted  for  different  values  of  the  time  t  differ,  therefore, 
only  by  a  factor  nearly  independent  of  height.  For  the  integrated 
profile  the  inflow  layer  is  less  well  defined  in  depth  and 
changes  in  the  course  of  one  period  in  the  interval  from  n  =  1 
to  n  s  2.  For  the  velocity  component  along  the  coast  the  ampli¬ 
tude  and  phase  functions  can  likewise  be  found.  The  hodograph 
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0.0  0.1  0.2  0.3 


0.0  0.1 

Fig.  2.7.  The  parameters  of  the  elliptical  wind  hodograph  as 
function  of  height  for  different  values  of  distance  from  the 
coastline.  The  full  lines  give  half  the  major  and  minor  axes, 
(-*-*-):  tQ,  ( - ):  tq  (fig.  2.6). 

The  amplitudes  refer  to  the  lower  scale,  and  the  phase 
angles  ♦0,t0  refer  to  the  upper  scale. 
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Fig.  2.7.  continued 


of  the  horizontal  velocity  vector  at  different  heights  above 
the  coast  are  given  by  |u( n) !cos( t+$u( n) ) ,  |v( n) |cos( t+$v( n) ) ) 
where  the  amplitude  and  phase  functions  corresponds  to  the  func¬ 
tions  for  two  components.  The  hodograph  is  an  ellipse  with  par¬ 
ameters  depending  on  height  and  distance  from  the  coast,  and  the 
behaviour  of  the  horizontal  wind  vector  can  therefore  be  charac¬ 
terized  by  the  parameters  of  this  ellipse  as  shown  in  Fig.  2.6, 
where  vertical  profiles  of  the  maximum  and  minimum  magnitude  of 
the  wind  vector,  together  with  the  phase  lag  of  the  maximum  mag¬ 
nitude  relative  to  the  maximum  land  surface  temperature,  and  the 
angle  between  the  maximum  wind  vector  and  the  gradient  of  sur- 
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face  temperature  are  plotted  for  different  distances  from  the 
coast . 


2.6.  Effects  of  horizontal  advection 


The  above  results  were  obtained  by  assuming  the  flow  to  be 
forced  solely  by  a  temperature  perturbation  at  the  surface.  A 
crude  estimate  of  the  effect  of  an  overlying  advective  velocity 
can  be  obtained  assuming  the  basic  state  velocity  to  be  constant. 
This  also  requires  that  it  is  permissible  to  neglect  the  wind- 
shear  near  the  surface  in  the  basic  state  velocity;  furthermore, 
the  neglect  of  the  "self  advection"  by  the  perturbation  u  in  com¬ 
parison  with  the  advection  by  the  background  velocity  U  means 
u  <<  U.  The  solution  in  this  case  can  be  obtained  directly  as 
above  by  noting  that  the  same  system  of  equations  and  boundary 
conditions  are  recovered  in  the  coordinate  system  moving  along 
with  the  velocity  U  in  the  x-direction.  To  be  more  specific,  it 
is  easily  seen  by  performing  the  Galilean  transformations,  that 
the  solution  in  the  original  (non-moving)  system  can  be  written 
for  each  wavenumber  in  real  notation  as; 


Re{u(k,n)  cos  k  ^  exp(-ix)}  5  u^  = 
k 

1  ( k)  ( k ) 

~2 1  a3+u)  <Z/H<S+J  cos(kS+T+*s+w(z/H3+(i)) ) 

+  AS-itz/Hu-J  cos(ke-T+i~lu(z/Hs_u))  } 


(2.35) 


where  A^^fn)  is  the  amplitude  of  the  u-component  corresponding 
to  the  wavenumber  k,  and  the  forcing  frequency  id  in  the  case  of 
U=0 .  Similarly,  ♦^Mn)  is  the  corresponding  phase  function, 

Hu  is  the  depth  scale  corresponding  to  the  frequency  «  ( Hw  ■ 

/K/w  ,  and  5  is  the  "induced"  forcing  frequency  ■  k  U/L.  The 
results  from  the  preceding  sections  are  recovered  for  3  «  0. 
Equation  (2.35)  enables  us  to  identify  two  limiting  cases  corre¬ 
sponding  to  small  and  large  scales  considered  above.  Assuming 
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a  >>  oj  corresponding  to  strong  advection  or  k  U/L  >>  w, 
Eq.  (2.35)  becomes 


(k)  (k) 

u(k)  =  A~  (z/Hj)  cos{ k£+$2  (2/H£))  cos( x) 


(2.36) 


The  scaling  introduced  in  the  preceding  sections  shows  that 
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On  introducing  2>  into  Eq.  (2.37)  we  obtain 
(k) 
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-here  ky  =  N/U  /k/w  .  Similarly  the  phase  function  can  be 
written  as 


(2.39) 


The  assumption  2  >>  i»  in  connection  with  diurnal  forcing,  where 
fg  ~  1  means  that  2  >>  f,  and  rotation  can  be  neglected.  The 
wavenumber  ky  defines  a  lower  limit  above  which  advection  by  U 
dominates  and  a!  >>  u>  is  equivalent  to  k  >>  ky.  With  this  assump¬ 
tion  the  first  argument  in  the  amplitude  function  in  Eq.  (2.38) 
is  small  compared  to  one  and  the  value  of  A  can  be  found  by  solv¬ 
ing  for  small  wavenumbers  as  in  the  previous  section  but  for 


fg  =  0. 
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For  the  phase  function  we  obtain 
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Using  these  two  expressions  valid  for  k  >>  ky  in  Eq.  (2.36)  the 
result  can  be  written: 


/ 


-  30  - 


=  i  exp[_ntf2^]  sin[i  kun]  sin[-k5+fc 


(2.42) 


Comparing  with  the  expression  obtained  for  large  wavenumbers  in 
the  case  of  zero  advection,  Eq.  (2.24),  the  depth  scale  is  de¬ 
creased  in  the  presence  of  advection  by  essentially  k1/6  k^1/2  « 
i^U.  For  large  scales  2  <<  ai  and  advection  have  only  little 
effect,  and  noting  that  4>u  =-$_  w  and  Atl)=A_  u  we  can  write 
Eq.  (2.35)  as: 


we  can  write 


/  u  \  ( k )  ( 

u'  =  A  (n)  cos (  kO  cos(t+i). 


(2.43) 


for  k  <<  ky.  This  is  identical  to  the  expression  valid  in  the 
case  u  =  0  as  expected.  For  intermediate  scales  k  »  ky  no  simple 
relation  exists  between  the  two  terms  on  the  right  in  Eq.  (2.35) 
and  each  wavenumber  must  be  treated  separately.  A  special  case 
appears  as  k  =  kD  equivalent  to  w=cd  ,  where  the  last  term  be¬ 
comes  undefined.  This  "resonance”  case  corresponds  (in  the  case 
of  zero  advection)  to  steady  convection  where  the  diffusive  depth 
scale  becomes  infinite.  The  proper  depth  scale  will  be  the  Ekman 
depth  =  /K/f  in  this  case.  Replacing  »  by  f  in  the  scaling  re¬ 
lations,  Eq.  (2.5),  the  set  of  equations  corresponding  to 
Eq.(2.13)  becomes 

32U  ,  i  , 
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-  u  =0 


w  =0 
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i  t  3w 
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The  equation  for  the  perturbations  pressure  p  becomes: 


(2.45) 


(The  prime  on  all  quantities  indicates  the  deviation  in  scaling). 
As  above,  the  solution  for  the  flow  can  be  obtained  from  Eq. 
(2.44)  and  the  boundary  conditions  with  the  exception  that  the 
boundary  conditions  at  infinity  for  v,  9,  and  p  cannot  be  satis¬ 
fied  because  of  the  appearance  of  a  linear  solution  in  Eq.  (244). 
Physically  this  particular  solution  corresponds  to  a  temperature 
perturbation  independent  of  height  with  associated  linear  press¬ 
ure  variation  and  linear  geostrophic  wind,  which  because  of  the 
orientation  of  the  coordinate  system  appears  only  in  the  v- 
component.  The  depth  of  the  cross  isobar  flow  u  is  determined  by 
the  Ekman  depth. 

The  effect  of  a  basic  state  velocity  will  be  pronounced  if  k^ 
is  not  large  compared  to  one  (U  ~  N  /k/  w  ):  Taking  reasonable 
values  for  the  parameters  viz.  N  =  0.01  s  ,  K  =  1  ms  , 
id  =  0.7  •  10~4  s-1  (=  1  day-1  )  then  the  requirement  is  CJ  ^ 

1  ms-1 .  The  smaller  scales  are  attenuated  at  all  levels  and 
phase  shifted  in  the  direction  determined  by  the  sign  of  U.  In 
the  above  arguments  we  have  assumed  U  to  be  positive;  negative 
U  yields  by  symmetry  identical  results  exept  that  -5  should  be 
replaced  for  £  in  Eq.  (2.35). 

It  is  well  known  that  the  inland  penetration  of  the  sea-breeze 
is  lessened  or  prevented  when  a  synoptic  scale  off-shore  wind 
is  present.  Empirical  predictions  of  the  occurrence  of  sea- 
breezes  have  been  based  on  the  negative  correlation  of  the  syn- 
optically  driven  wind  and  sea-breeze  occurrence  (Biggs  and 
Graves,  1962).  Based  on  a  linear  model  Walsh  (1974)  finds  that 
this  model  lends  theoretical  support  to  the  empirical  relations. 
His  analysis  is  based  on  the  assumption  that  the  requirement  for 
the  occurrence  of  sea-breeze  is  umax  2  |U|  where  umax  is  the 
maximum  velocity  predicted  by  the  linear  model.  Clearly  this  re¬ 
quirement  is  inconsistent  with  the  linearization  of  the  advection 
terms.  In  addition,  large  amplitudes  of  the  surface  heating 
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necessary  to  meet  the  above  requirement  for  realistic  values  of 
the  windspeed  U  and  the  background  stratification  N2  leads  to 
violation  of  the  linearization  of  the  adiabatic  heating  term  in 
the  thermodynamic  equation.  The  necessary  assumption  here  is 
that  the  flow-induced  perturbation  of  the  stratification  can  be 
neglected.  This  means  that  the  following  relation  must  hold: 


A  9 


AN' 


g 

H  0 


(2.46) 


.he  analysis  by  Walsh  (1974)  is  based  on  a  particular  choice  for 
N  since  he  did  not  anticipate  the  possibility  of  scaling  the 
problem  such  that  only  fs  =  f/w  appears  as  an  external  parameter 
as  shown  in  the  present  analysis.  Using  his  results  of  the 
maximum  value  umax  in  the  linear  model  as  function  of  a  back¬ 
ground  velocity  it  is  possible  to  obtain  umax  as  function  of 
the  wavenumber  k0  as  shown  in  Fig.  2.7.  The  results  are  well 
reproduced  by  the  relation 
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(2.47) 


with  a  *  0.65.  Eq.  (2.47)  reproduces  the  asymptotic  behaviour  for 
umax  as  U  +  0  ( ky  ♦  <*>)  and  also  as  D  ♦  *  (ky  ♦  0):  umax  ♦  0  as 
1/U.  At  moderate  values  of  ky  Fig.  2.7  shows  Eq.  (2.47)  to  re¬ 
produce  the  model  results  adequately. 


Returning  to  the  requirement  Eq.  (2.46)  in  connection  with  the 
assumption  umax  *  U  (where  denotes  dimensional  quantity): 


g  A0o  (U)  g  46o  (0) 
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for  kn  <<  1  we  obtain 
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since  our  linearization  requires  /N  <<  1  and  umax  a  =  0.14 
clearly  umax  *  U  strongly  violates  the  linearization  for  U 
large.  Ignoring  for  a  moment  this  objection  Eq.  (2.49)  can  be 
rewritten  as 

o  \rr  o 

A 0  >  -  U*  (2.50) 

0  “  0.14  g  Ik 

This  relation  is  similar  to  the  expression  found  by  Walsh  (1974) 
and  has  the  form  of  the  empirical  relation  found  by  Biggs  and 
Graves  (1962);  furthermore,  as  shown  by  Walsh  the  consta  t  in 
Eq.  (2.50)  is  very  nearly  the  empirical  constant  when  a  realistic 
value  for  K  is  inserted* 

In  view  of  the  simplicity  of  the  linear  model  and  the  objections 
raised  above  it  is  surprising  that  the  model  agrees  with  obser¬ 
vational  data  with  an  accuracy  as  good  as  that  with  which  the 
parameters  in  the  model  can  be  determined.  The  appropriate  mean 
diffusivity  K  cannot  be  estimated  independently  better  than  as 
an  order  of  magnitude  estimate;  also  the  relation  between  the 
gradient  wind,  which  enters  in  the  empirical  relation  similar  to 
Eq.  (2.49)  in  place  of  U,  and  the  "bulk"  background  velocity  U 
used  in  the  model  must  include  a  geometrical  factor  of  order  one 
to  account  for  the  effect  of  windshear  in  the  boundary  layer. 

In  the  case  of  zero  background  velocity  the  linear  model  predicts 
the  magnitude  of  the  maximum  cross  isobar  flow  umax  to  be  pro¬ 
portional  to  the  amplitude  of  the  land  surface  temperature.  The 
scaled  profile  of  u  obtained  from  the  linear  model  by  solving 
for  all  wavenumbers  and  integrating  (Eq.  2.8)  is  shown  in  Fig. 

2.8  for  different  values  of  t  =  u  t.  The  profile  corresponding 


*)  The  numerical  constant  (  *  0.14)  in  Eq.  (2.49)  differs  from 
the  constant  found  by  Walsh  (0.11)  because  of  the  fitting  used 
here  by  the  expression  Eq.  (2.46).  Using  (2.49)  and  the  empirical 
coefficient  yields  K  ■  8  m^s”^ . 
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Fig.  2.8.  The  effect  on  the  maximum  value  of  the  u-component 
of  an  background  velocity  U.  ky=N/U/K/u>.  (o):  Results  from 
rescaling  data  from  Walsh(1974),  ( - ):  Eg. (2. 47). 


to  t  =  0.15  corresponds  to  the  profile  when  the  largest  value 
of  u  is  found.  The  model  estimate  of  umax  is  therefore 


~(0) 

umax 


0.22 


(2.51  ) 


where  the  superscript  zero  indicates  zero  advection.  The  numeri¬ 
cal  constant  is  valid  only  for  the  case  of  a  step  change  in  tem¬ 
perature  across  a  straight  coastline  which  is  considered  here, 
and  also  only  valid  for  the  particular  choice  of  the  nondimen- 
sional  Coriolis  parameter  ( f s  =  f/w  =  1.5)  corresponding  to 
diurnal  forcing  at  mid  latitudes.  The  strong  dependence  found 
above  on  the  background  velocity  and  the  favorable  comparison 
with  empirical  data  suggest  that  a  reasonable  order  of  magni¬ 
tude  estimate  of  the  maximum  velocity  taking  into  account  the 
self-advection  at  finite  values  of  A60  can  be  obtained  using 
Eq.  (2.47),  and  assuming  the  effective  advective  component  to  be 
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proportional  to  umax ,  viz  U  =  a  umax»  where  a  is  a  geometrical 
factor  of  order  one.  The  proposed  equation  becomes 


umax 


(0) 

umax 


( l-exp(-aky) ) 


(2.52) 


..here  the  last  equality  follows  from  the  definition  of  ky.  Eq. 
(2.51)  can  be  solved  for  ky  by  iteration  and  umax  is  then  de¬ 
termined  by  the  last  expression.  From  Eq.  (2.50)  it  is  clear 
that  the  linear  result  becomes  questionable  as  60o  increases  and 
as  stratification  N  decreases.  Using  Eq.  (2.52)  in  the  limit 
ky  <<  1  we  obtain 


umax  ■ 

where  umax 
yields 

“max 


t  9Uma* 


=  0.22,  a  =0.65.  Using  a  =  a  =0.65  and  K 
24  «  (— ) 


(2.53) 


2  -1 
m  s 


(2.54) 


Fig.  2.9.  ( - ):  from  the  self  advection  model  from  the 

solution  of  Eq.  (2.52)  with  a-a*0.65.  ( - ):  Linear  model. 
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Comparing  Eq.  (2.53)  with  (2.51)  we  see  that  the  strong  depen¬ 
dence  on  stratification  in  the  linear  model  disappears  and  is  re¬ 
placed  by  a  very  weak  dependence  on  the  diffusivity  K  in  (2.53). 
Equation  (2.53)  can  be  interpreted  as  a  Froude  number  relation¬ 
ship: 

/  Ap\_1/2 

=  “max  {<3*  -1.5  (2.55) 

where  d  is  the  depthscale  =  /K/  uT  ,  ip/p0  =  &6o/0,  and  we  have 
used  a  =  a  =  0.65.  In  studies  of  the  sea-breeze  front  Simpson  et 
al.  (1977)  have  found  that  if  d  is  the  depth  of  the  intruding 
(cold)  sea  air,  u__v  is  the  rate  of  advance  of  the  sea-breeze 

ilia  A 

front  and  Ap  is  the  density  difference  between  sea-  and  land  air 
then  Fr  s  0.7.  In  our  definition  u__„  represents  the  maximum 
air-speed  and  comparing  the  two  Froude  numbers  we  find  that  they 
agree  if  the  rate  of  advance  of  the  sea-breeze  front  is  2.1 
times  slower  than  the  maximum  speed.  This  number  agrees  sur¬ 
prisingly  well  with  what  is  found  from  observations.  As  an  ex¬ 
ample  Kimura  and  Eguchi  (1978)  find  the  average  ratio  to  be  2.2 
based  on  an  average  over  56  sea-breeze  days  observed  from  13 
surface  stations  in  Japan.  The  maximum  speed  predicted  from 
Eq.  (2.54)  can  be  likewise  compared  with  observations;  in  the 
field  study  of  the  sea-breeze  at  the  great  lakes  in  the  United 
States  Keen  and  Lyons  (1978)  found  a  maximum  onshore  windspeed 
of  4.1  ms-1  at  a  height  of  150  m  above  the  shoreline,  the  tem¬ 
perature  difference  between  land  and  sea  was  8  K.  Using  this 
value  Eq.  (2.54)  yields  umax  =  4.0  ms-^  .  Compare  this,  however, 
with  the  variability  discussed  in  connection  with  the  numerical 
simulations  (Chapter  3). 


2.7.  Comparison  with  inviscid  model 

The  result  from  the  linear  model  for  the  cross  isobar  flow  has 
already  been  discussed  in  the  preceding  section.  The  v-com- 
ponent  at  the  coastline  at  low  heights  is  much  smaller  than  the 
u-component  in  magnitude,  reflecting  the  influence  of  high  wave- 
numbers  where  rotation  is  negligible.  The  hodograph  of  the  hori- 
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zontal  wind  is  an  ellipse  in  all  points  because  of  the  assumed 
periodicity.  The  assymetry  observed  for  the  land  and  sea  breezes 
can  partly  be  incorporated  into  the  linear  model  by  adding  har¬ 
monic  terms  to  simulate  the  typical  diurnal  land  surface  tem¬ 
perature.  As  discussed  by  Mak  and  Walsh  (1976)  this  will  account 
for  only  part  of  the  assymmetry  because  a  large  part  is  caused 
by  the  diurnal  variation  of  stratification  and  diffusivity  which 
cannot  be  incorporated  into  the  linear  model.  Because  of  the 
strong  influence  of  the  small  scales  on  the  magnitude  of  the 
maximum  cross  isobar  component  this  quantity  is  to  only  a  rela¬ 
tively  small  degree  influenced  by  complicated  initial  conditions 
at  the  beginning  of  a  sea-breeze  day.  A  slowly  varying  synoptic 
scale  u-component  will  act  to  modify  the  maximum  breeze  u- 
component  as  described  above.  The  v-component,  and  also  the  u- 
component  at  higher  levels,  is  dominated  by  larger  scales  where 
the  Coriolis  force  is  important  and  the  initial  conditions  will 
likewise  be  important.  At  all  points  the  wind  vector  exibits  the 
expected  clockwise  turning  (for  f  >  0). 

Near  the  coastline  very  large  gradients  exist  in  the  vertical 
velocity,  which  is  antisymmetric  across  the  coastline.  The  be¬ 
haviour  of  w  here  is  determined  by  the  contribution  from  large 
wavenumbers  for  which  advection  becomes  important  even  for  small 
values  of  the  forcing.  As  discussed  in  the  previous  section  the 
u-field  also  is  dominated  by  large  wavenumbers  near  the  coastline. 
Consequently,  the  detailed  results  from  the  linear  model  cannot 
be  expected  to  apply  to  the  land-  and  sea-breeze  near  the  sea- 
breeze  front.  If  we  define  the  "inner  region"  in  the  circulation 
as  the  region  where  the  non  linear  effects  are  important,  then 
we  can  expect  the  linear  dynamics  to  apply  in  the  complementary 
"outer"  region.  The  reason  that  such  a  division  is  meaningful  is 
that  we  are  dealing  with  a  disturbance  generated  essentially  at 
one  point  and  propagating  in  a  dispersive  medium.  The  propagation 
of  a  disturbance  in  a  stratified  inviscid  fluid  is  usually  ana¬ 
lysed  in  terms  of  internal  gravity  waves,  and  such  an  analysis 
applied  to  the  sea  breeze  outer  region  was  presented  by  Geisler 
and  Bretherton  (1968)  who  treated  the  problem  as  an  initial-value 
problem  by  introducing  the  temperature  distribution  everywhere 
instantaneously.  They  introduced  the  term  the  sea-breeze  fore- 
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runner  as  the  disturbance  being  carried  by  the  gravity  waves  and 
reaching  points  away  from  the  coastline  in  advance  of  the  sea- 
breeze  proper.  In  their  analysis  rotation  and  diffusion  were 
neglected,  and  they  found  in  the  simplest  case  where  stratifi¬ 
cation  is  independent  of  height  that  the  perturbation  in  Q  can 
be  expressed  as  g/N(z/h,T)  where  h  is  the  depth  of  the  initially 
imposed  temperature  perturbation  and  t  =  Nh  t/x.  The  connection 
to  the  linear  harmonic  model  can  be  seen  if  h  is  interpreted  as 
the  diffusive  depth  scale  /k/uj  .  The  dispersion  relation  for 
gravity  waves  in  the  system  is  given  by  (Pearson,  1973): 

i»2  =  f2  +  N2k2/m2  (2.56) 

when  rotation  is  included.  Neglecting  rotation  the  horizontal 
phase  and  group  speed  for  waves  with  vertical  wavenumber  m  is 
u>/k  =  N/m,  and  noting  that  all  modes  are  excited  simultaneously 
as  the  temperature  contrast  is  imposed  in  the  Geisler  and 
Bretherton  model,  at  any  point  away  from  the  coastline  a  sequence 
of  modes  arrive  with  the  largest  vertical  scales  arriving  first 
and  eventually  modes  of  smaller  and  smaller  vertical  wavelength 
arrive.  In  the  periodic  model  three  other  mechanisms  are  in  ef¬ 
fect,  firstly  rotation  is  included,  affecting  the  behaviour  of 
the  large  scales;  secondly,  the  periodic  forcing  means  that  the 
disturbance  arriving  at  a  point  is  a  superposition  of  disturb¬ 
ances  generated  at  different  times  with  a  strength  given  by  the 
harmonic  time  variation;  and  thirdly,  the  effect  of  surface  drag 
and  the  internal  dissipation  is  to  damp  the  different  modes  at 
different  rates.  The  e-folding  time  on  modes  with  vertical  wave- 
number  m  is  of  order  m-2  K"1  which  for  modes  of  vertical  extent 
of  h  ■  /K/w  m-1  becomes  equal  to  L  =  N/w/K/io,  which  is  the 
horizontal  scale  length  in  the  periodic  model.  The  largest  per¬ 
turbation  in  u  occurs  at  the  surface  in  the  inviscid  model  where 
a  slip-condition  is  assumed.  Employing  the  scaling  introduced  for 
the  periodic  model,  the  result  found  by  Geisler  and  Bretherton 
for  u(x,z,t)  can  be  written  as 
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(2.57) 


where  we  have  used  Eq.  (24)  in  the  original  paper  and  h  =  /k/uj, 
L  =  N/ui/k/w,  etc.  from  the  scaling  used  in  the  periodic  model 
described  above.  It  may  seem  confusing  to  introduce  a  forcing 
frequency  u>  and  a  diffusivity  K  when  we  are  dealing  with  an 
inviscid  model  with  steplike  forcing,  however  it  makes  compari¬ 
son  with  the  periodic  model  easier,  as  will  be  apparent  below. 
At  the  surface  the  result  Eq.  (2.57)  can  be  written: 
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(2.58) 


This  function  is  shown  in  Fig.  2.10.  It  can  be  seen  that  the 
maximum  value  found  in  the  periodic  model  (20. 22)  is  exceeded 
already  at  values  of  r/Z  ~  2,  suggesting  the  importance  of 
diffusion  of  momentum.  Introducing  into  the  model  an  internal 
dissipation  of  the  waves  leads  to  the  appearance  of  the  term 
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Fig.  2,10.  The  u-component  at  the  surface  from  the  inviscid 
model  of  Geisler  and  Bretherton  (Eq.  2.58). 
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5=x/L=  x/(h3N/K) 

Fig.  2.11.  The  u-component  tn  the  Geisler  and  Bretherton  model, 
in  the  case  where  internal  dissipation  is  included.  Numbers 
on  the  curves  refer  to  the  time  x. 

exp(-y3  5)  inside  the  integral  in  Eq.  (2.57),  as  also  discussed 
in  the  paper  by  Geisler  and  Bretherton.  The  effect  of  this  term 
is  to  give  a  sharp  cutoff  at  y  =  £-1/3.  The  maximum  value  of 
u  now  becomes  approximately  equal  to  Eq.  (2.58)  with 
replaced  for  x/5.  The  result  is  shown  in  Fig.  2.11.  The  curve 
labelled  x  =  »  gives  the  asymptotic  value  of  umax,  and  the  curves 
labelled  with  different  valves  of  x  gives  u„,ax  (5,x,0).  These 
curves  follow  the  asymptotic  curve  up  to  the  point  where  the 
diffusive  cutoff  occurs  at  5  =  x3/2.  Note  that  we  can  eliminate 
the  dependence  on  w  by  substituting  w  =  K/h2  with  K  -  x/(Nh3/K) 
and  x  =  t/(h2/K).  The  maximum  u-component  in  the  periodic  model 
is  shown  together  with  the  x  =  »  curve  from  Fig.  2.11  in  Fig. 

2.12.  The  large  difference  is  a  consequence  of  the  mechanisms 
active  in  the  periodically  forced  model:  the  surface  drag  which 
becomes  infinite  at  5  *  0  (cf.  Eq.  2.24),  the  finite  tempera¬ 
ture  gradient  above  the  surface  at  4  *  0  resulting  from  the 
action  of  the  vertical  velocity  in  the  thermodynamic  equation, 
and  the  periodic  forcing.  The  horizontal  temperature  distri¬ 
bution  in  the  periodic  model  is  illustrated  in  Fig.  2.13  for 


5=x/L 

Fig .  2.12.  The  maximum  u-component  as  function  of  distance 
from  the  coast  for  the  periodic  model  (lower  curve),  and 
the  Geisler  and  Bretherton  model  with  internal  dissi¬ 
pation  (curve  marked  t  =  °°  on  the  previous  figure). 


Fig .  2.13.  The  temperature  perturbation  at  the  time  of  maximum 
perturbation  at  the  surface  in  the  periodic  model  for  differ¬ 
ent  values  of  height.  (The  zero  point  for  8  has  been  chosen 
arbitrarily  for  each  curve.) 
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different  values  of  n  at  the  time  of  maximum  in  the  surface 
temperature  difference.  The  curves  demonstrate  the  difference 
in  the  forcing  function  between  the  viscous  periodic  model  and 
the  model  by  Geisler  and  Bretherton  where  the  initial  forcing 
is  a  stepfunction  of  amplitude  equal  to  one  up  to  the  height 
n  =  1  and  zero  above.  It  should  be  noted  that  the  vertically 
integrated  temperature  perturbation  far  away  from  the  coast  is 
identical  in  the  two  models.  The  effect  of  surface  drag  can  be 
incorporated  into  the  model  of  Geisler  and  Bretherton,  but  only 
in  a  crude  way,  as  discussed  in  their  paper.  In  the  periodic 
model  treated  sofar  the  atmosphere  is  modelled  as  one  layer  of 
constant  diffusivity  and  basic  state  stratification.  This  is 
obviously  a  gross  oversimplification  and  even  in  the  outer  part 
of  the  flow  large  deviations  in  the  structure  must  be  expected, 
in  general,  because  of  temporal  and  spatial  structure  in  the 
basic  state  parameters.  The  most  obvious  inhomogeneity  is  be¬ 
tween  the  boundary  layers  developed  over  sea  and  land  because 
of  the  difference  in  surface  heating.  In  the  outer  part  one  may 
expect  the  conditions  over  land  and  sea  to  be  uncoupled  and  the 
inviscid  model  by  Geisler  and  Bretherton  (and  the  periodically 
forced  inviscid  model  described  by  Kimura  and  Eguchi  (1978)) 
could  apply  in  particular  over  the  sea,  whereas  for  conditions 
over  land  turbulent  diffusion  must  be  incorporated.  Over  land 
there  may  be  no  outer  region  in  the  boundary  layer  but  only 
above  it  because  as  one  moves  inland  from  the  coast  in  the  con¬ 
vective  boundary  layer  the  stratif ication  vanishes  except  near 
the  surface  where  it  is  superadiabatic  in  the  sea-breeze  situ¬ 
ation.  As  a  consequence,  the  phase  speed  of  gravity  waves  tend 
to  zero  and  the  linear  dynamics  breaks  down.  In  the  periodic 
viscous  model  the  determination  of  the  proper  diffusivity  presents 
a  problem,  since  it  has  a  large  dependency  on  height.  The  values 
used  in  the  example  above  are  proper  only  for  near  surface 
conditions.  The  stratification  enters  only  in  the  thermodynamic 
equation  through  the  adiabatic  cooling  term  w  •  r  which  is  un¬ 
important  near  the  surface  but  is  responsible  for  the  cooling 
over  land/heating  over  the  sea  at  heights  in  the  upper  part  of 
the  boundary  layer  over  land  and  above,  generating  high  pressure 
over  land  and  the  return  flow.  The  stratification  parameter  N  in 
the  linear  model  should  therefore  properly  be  a  characteristic 
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value  for  the  upper  boundary  layer  including  the  (possible) 
inversion  and  the  stratification  above.  With  an  inversion  on  top 
of  a  mixed  layer  both  the  turbulent  diffusivity  and  the  strati¬ 
fication  depends  on  height,  and  the  velocity  and  temperature 
profiles  will  deviate  from  the  profiles  found  in  the  simple 
model.  The  effect  of  an  elevated  inversion  layer  on  the  profiles 
in  the  framework  of  the  linearized  periodically  forced  model  is 
discussed  in  the  next  section. 


2.8.  Vertical  variation  of  diffusivity  and  stratification 

In  real  boundary  layers  the  diffusivity  of  heat  and  momentum 
are  typically  functions  of  height  through  the  dependence  of  the 
mean  flow  structure  and  the  stratification.  The  inclusion  of  a 
height-dependent  diffusivity  and  stratification  in  the  linearized 
model  leads  to  an  eight-order  differential  equation  with  non¬ 
constant  coef f icients,  and  no  simple  solutions  exist  in  general. 
Semi  analytical  solutions  can  be  obtained,  however,  if  the  par¬ 
ameters  can  be  assumed  to  be  piecewise  constant.  For  each  layer 
in  which  K  and  N  are  constant  the  solution  is  identical  to  the 
solutions  considered  above  except  for  the  boundary  conditions. 

The  total  solution  is  obtained  by  applying  matching  conditions 
at  the  interfaces  between  layers  in  addition  to  the  boundary 
conditions  at  the  surface  and  at  infinity.  The  matching  conditions 
are  that  all  fields  (u,  v,  w,  B  and  £>)  and  the  fluxes  of  heat 
and  momentum  (K  3u/3z,  K  3v/3z,  K  3$  3z)  are  continuous  across  the 
interfaces.  Subdividing  into  N  homogeneous  layers  gives  ( N— 1 ) 
interfaces,  and  with  our  previous  boundary  conditions  at  the  sur¬ 
face  (u*v=w=0  and  0  prescribed)  together  with  the  boundary  con¬ 
ditions  at  infinity  (four  coefficients  to  exponentially  increas¬ 
ing  terms  vanish)  we  have  8  N  constraints  to  determine  the  8  N 
coefficients  in  the  pressure  profile.  The  other  fields  are  found 
from  the  pressure  profile  by  applying  the  relations  (2.22)  for 
each  layer.  The  large  number  of  equations  appearing  when  we  sub¬ 
divide  prohibits  any  manageable  explicit  expressions  for  the  ver¬ 
tical  profiles,  but  since  the  coefficient  matrix  is  sparse  the 
solution  is  easily  obtained  numerically  even  in  the  case  of  many 
sublayers. 
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With  discontinuties  in  the  diffusivity  K  this  means  that  the 
fields  will  not  have  a  continuous  first  derivative  across  the 
interfaces;  note,  however,  that  the  first  derivative  of  the  ver¬ 
tical  velocity  and  the  pressure  are  continuous  by  virtue  of  the 
matching  conditions  together  with  the  continuity  equation  and 
the  hydrostatic  equation,  respectively.  The  matching  conditions 
ensure  non-singular  behaviour  of  all  terms  in  the  dynamical 
equations  when  applied  across  interfaces. 

With  this  model  we  can  simulate  the  effect  of  an  elevated  inver¬ 
sion  layer  by  subdividing  into  3  layers.  If  the  height  of  the 
inversion  base  is  larger  than  the  diffusive  scale  height  in  the 
lower  layer  (*/K;/(d)  we  must  expect  the  scaling  Eqs.  (2.5)  to 
be  applicable  but  with  the  flow  modified  by  the  inversion.  If  the 
inversion  is  at  a  height  lower  than  this  scale  height  then  no 
unique  scaling  exists  and  the  structure  and  strength  of  the  flow 
is  a  function  of  the  basic  state  parameters  in  all  three  layers. 
If,  on  the  other  hand,  the  inversion  base  is  at  a  height  much 
larger  than  the  scale  height  the  effect  of  the  inversion  is  neg- 
negligible. 

For  scaling  the  additional  non-dimensional  parameters  the  par¬ 
ameters  in  the  lowest  layer  are  used.  This  the  height  of  the 
inversion  base,  inversion  thickness,  the  scale  heights  and  the 
scale  lengths  in  the  two  upper  layers.  The  proper  diffusivity 
for  each  layer  is  strongly  dependent  on  the  stratification  and 
a  convective  boundary  layer  is  generally  capped  by  an  inversion 
layer  effectively  reducing  diffusive  transport  above  the  inver¬ 
sion.  With  horizontal  homogeneous  conditions  the  inversion  layer 
is  formed  by  the  entrainment  process  at  the  top  of  the  boundary 
layer.  The  strength  will  depend  on  the  stratification  above  the 
inversion,  the  rate  of  surface  heating,  and  possible  subsidence 
on  the  synoptic  scale.  In  the  sea-breeze  situation  conditions 
are  complicated  by  the  inhomogeneity  across  the  coastline,  the 
effect  on  stratification  by  the  vertical  motions  generated  by 
the  sea-breeze  itself  and  horizontal  advection.  These  compli¬ 
cations  cannot  adequately  be  incorporated  into  the  linear  model, 
and  presumably  they  can  be  handled  only  numerically.  In  spite  of 
this  it  is  of  interest  to  study  the  effect  of  an  inversion  layer 
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on  the  flow  in  the  linear  model  since  the  linear  dynamics  still 
apply  in  the  circulation  locally  with  the  possible  exception 
of  the  region  near  the  sea-breeze  front  and  the  coastline.  With 
essentially  vanishing  diffusivity  above  the  inversion  the  depth 
scale  of  the  flow  becomes  very  small  and  the  flow  in  the  boundary 
layer  will  be  influenced  only  by  the  thermal  stratification  up 
to  a  height  immidiately  above  the  inversion  base.  This  means  that 
the  effect  of  an  inversion  in  the  framework  of  the  linearized 
model  can  be  modelled  using  only  two  layers. 

The  maximum  magnitude  of  the  two-level  wind  vector  is  shown  as 
function  of  the  ratio  of  the  Brunt-Vaisala  frequencies  for  dif¬ 
ferent  heights  of  the  inversion  base  on  Figs.  2.14a-c.  The  depth 
scale  in  the  upper  layer  is  set  to  0.1  times  the  depth  scale  in 


N2/N, 


Fig.  2.14.  Maximum  value  of  the  low  level  speed  as  function 
of  the  ratio  of  the  values  for  N  in  the  two  layers. 

Values  given  on  the  curves  refer  to  the  height  of  the  lower 
layer  (inversion  height)  scaled  with  parameters  pertaining 
to  the  lower  layer.  K  on  the  figures  gives  the  wavenumber  k. 
(See  text) . 
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2.14.  continued 


the  lower  layer  corresponding  to  a  decrease  of  the  diffusivity 
by  a  factor  of  100  across  the  inversion  base.  This  value  is  some¬ 
what  arbitrary,  but  using  larger  values  of  this  ratio  results  in 
only  very  small  changes  in  the  maximum  low-level  wind  speed.  As 
expected  from  the  discussion  above,  however,  the  depth  of  the 
flow  above  the  inversion  base  is  reduced  as  the  depth  scale  is 
decreased.  For  all  wavenumbers  the  wind  speed  is  strongly  influ- 
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enced  by  the  reduced  vertical  mixing  even  in  the  absence  of  tem¬ 
perature  inversion.  The  values  for  the  homogeneous  case  (h  =  00  ) 
are  close  to  the  values  for  h  *  2  and  there  occurs  obviously  an 
"overshoot"  for  values  of  h  comparable  to  one.  For  larger  ratios 
of  the  stratification  the  wind  speed  decreases  toward  an  asymp¬ 
totic  value  determined  by  the  inversion  height,  the  flow  becomes 
confined  to  the  lower  layer  except  for  a  thin  layer  with  a  strong 
return  flow  in  the  inversion  layer.  The  phase  angle  ( tq  on  Fig. 
2.6)  decreases  as  the  inversion  strength  is  increased  and  with  a 
strong  low  level  inversion  t0  ~  0  for  all  wavenumbers,  the  angle 
<t>0  also  tends  to  zero  for  strong  inversions.  This  is  not  surpris¬ 
ing  in  view  of  the  discussion  in  section  2.7  and  the  results  in 
section  2.4;  the  stratification  is  most  important  in  the  dynamics 
near  n  =  1  (in  the  homogeneous  case),  whereas  the  largest  effect 
of  diffusion  is  near  the  surface;  consequently  a  two-layer  model 
with  an  inversion  layer  at  n  -  1  must  dynamically  become  similar 
to  the  one-layer  model  for  large  wavenumbers.  We  have  used  the 
parameters  in  the  lowest  layer  for  scaling,  and  if  the  "effec¬ 
tive"  stratification  increases  with  the  strength  of  the  inver¬ 
sion,  then  for  any  wavenumber  (=  where  i  is  the  physical 

wavelength  and  L-j  is  the  scale  length  in  the  lower  layer  *  ) 

the  "effective"  wavenumber  influenced  by  the  much  larger  length 
scale  L2  above  the  inversion  base  must  be  larger  than  the  wave- 
number  itself.  Figure  3.2  shows  windprofiles  for  the  more  compli¬ 
cated  situation  where  we  have  employed  several  layers  simulating 
a  parabolic  diffusivity  profile  and  a  jump  in  the  Brunt-Vaisala 
frequency  at  the  top  of  the  layer  of  strong  mixing.  The  increase 
of  the  diffusivity  with  height  near  the  ground  tends  to  concen¬ 
trate  a  strong  windshear  near  the  ground  to  lower  the  level  of 
maximum  wind  speed  compared  to  the  homogeneous  case,  the  return 
flow  is  concentrated  near  the  top  of  the  layer  with  strong 
mixing,  as  expected  from  the  discussion  above. 

In  connection  with  an  inversion  layer  with  an  effective  decoup¬ 
ling  of  the  horizontal  wind  a  synoptic  pressure  gradient  can 
give  rise  to  a  strong  wind  shear  across  the  inversion,  and  the 
role  of  horizontal  advection  by  such  a  basic  state  wind  velocity 
is  complicated  by  the  different  advection  rates  acting  on  the 


Fig.  2,15a.  Contour  plot  of  the  u-component  at  the  time 
of  maximum  surface  temperature  perturbation  as  function  of 
%  (abcissa)  and  n  (ordinate).  Contour  interval  is  0.025. 
(Max  value  is  7*0.025=0.18). 


Fig,  2.15b.  As  fig.  2.15a,  but  for  the  case  with  background 
velocity  0  with  kjj*4. 


Fig.  2.16a.  As  fig.  2.15a,  but  in  the  case  of  a  two  layer 
structure  (inversion)  with  N2  in  the  layer  above  n=1  given 
as  N2=100*Nlf  and  K2=0.01*K^.  Scaling  with  lower  layer 


Fig.  2.16b.  As  fig.  2.16a,  but  with  a  background  U  with  kg=4 
in  both  layers. 
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flow  in  the  mixed  layer  and  on  the  return  flow  near  the  inver¬ 
sion.  In  consistency  with  the  treatment  of  a  constant  horizontal 
basic  state  velocity  (section  2.6)  and  with  the  layered  model 
treated  here,  different  basic  state  velocities  can  be  assigned 
to  each  layer  to  simulate  this  situation.  For  the  two-layer  model 
including  windshear  at  the  inversion  results  are  shown  on  Fig. 
2.16  for  the  cross  isobar  flow  at  the  time  of  maximum  surface 
temperature . 

The  effect  of  the  inversion  in  concentrating  the  return  flow 
and  reducing  the  low  level  wind  speed  is  apparent  from  these 
figures;  the  effect  of  a  constant  basic  state  velocity  is  quali¬ 
tatively  the  same  with  and  without  the  inversion.  At  the  same 
distance  inland  the  windspeed  actually  increases  with  increasing 
basic  state  velocity  (decreasing  ky)  because  the  maximum  is 
displaced  away  from  the  coast  in  the  direction  of  the  wind.  The 
maximum  is  decreased,  however,  as  discussd  in  section  2.6. 

The  effect  of  different  basic  state  velocity  in  the  two  layers 
is  to  displace  the  points  of  maximum  low-level  inflow  and  maximum 
return  flow  in  the  inversion  layer  away  from  each  other  in  the 
horizontal.  With  larger  shears  the  situation  changes  and  disturb¬ 
ances  generated  in  the  shear  layer  appear  equivalent  to  the  un¬ 
stable  gravity  waves  which  are  generated  at  the  interface  be¬ 
tween  two  homogeneous  layers  of  different  density  and  with  a 
velocity  jump  across  the  interface  (see,  e.g.  Gossard  and  Hooke, 
1975).  The  stability  analysis  is  complicated  in  the  case  treated 
here  by  the  surface  boundary  conditions  and  the  diffusivity  in 
the  two  layers.  In  the  model  considered  where  all  fields  are 
assumed  harmonic  in  time  with  frequency  w,  and  where  9  is  tied 
at  the  surface,  these  waves  appear  as  finite  disturbances  which 
are  amplified  with  increasing  shear.  This  is  obviously  a  model 
artifact  and  their  structure  in  the  real  atmosphere  cannot  be 
analysed  in  the  framework  of  this  simple  model.  Waves  are  often 
observed  in  connection  with  the  nocturnal  boundary  layer  where 
they  may  form  intermittently  when  strong  shear  is  present.  In 
the  laboratory  experiments  by  Britter  and  Simpson  (1976)  such 
waves  are  found  behind  the  sea-breeze  head  on  top  of  the  layer 


52 


of  protruding  sea  air  simulated  by  a  gravity  current  of  a  dense 
fluid  undercutting  a  lighter  fluid. 


2.9.  Circular  island 


In  the  previous  sections  we  have  considered  the  two-dimensional 
problem  of  flow  across  a  straight  coastline.  In  cases  where  the 
coastline  is  curved  the  flow  will  be  subject  to  imposed  diver¬ 
gence  or  convergence  which  in  turn  will  induce  vertical  motions 
and  pressure  adjustments.  Because  of  the  1/N  dependence  of  hori¬ 
zontal  velocity  which  is  coupled  to  the  thermodynamic  equation 
through  the  adiabatic  tooling  term  w r,  the  expected  effect  of 
increased  convergence/divergence  is  to  decrease  the  magnitude 
of  horizontal  velocity  compared  to  the  case  of  the  straight 
coastline. 


The  effect  can  be  studied  by  solving  the  linear  flow  problem  for 
the  case  of  a  heated  circular  island.  Assuming  angular  symmetry, 
and  interpreting  and  horizontal  coordinate  5  as  the  (scaled) 
radial  distance  from  the  island  center,  the  dynamical  equations 
Eq.  (2.6)  are  unaltered  except  for  the  continuity  equation,  which 
in  cylindrical  coordinated  becomes: 


u  3u  3w 

_  +  _  +  _ 

C  3£  3  n 


0 


(2.59) 


The  Fourier  representation  is  impractical  in  this  case  and  the 
fields  are  instead  expanded  in  a  Fourier-Bessel  series: 


u(  S,ti)  *  -  l  u(l,n)  J 1  ( 0,5/D)  •  0  <  5  <  D  (2.60) 

for  the  radial  velocity  component  u  and  similarly  for  the  tan¬ 
gential  velocity  v.  The  temperature,  vertical  velocity,  and 
pressure  are  given  the  expansion,  viz  for  9: 

0» 

-  I  e(t,n)  Jo(0,5/D)  t  0  <  5  <  D 

t-1 


e(C,n) 


(2.61) 
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The  functions  JQ  and  are  Bessel  functions  of  zero'th  and 
first  order,  respectively,  and  aj  is  the  t'th  zero  of  JQ.  Inser¬ 
ting  the  expansions  (2.60)  and  (2.61)  into  Eqs.  (2.6)  with  the 
continuity  equation  replaced  by  Eq.  (2.59)  leads  as  before  to  Eqs 
(2.13)  for  the  expansion  coefficients.  The  wavenumber  k  now  takes 
the  discrete  values  kj  =  oa/d,  l  =  1,...,  ».  In  particular  we 
can  rewrite  the  continuity  equation  for  the  i'th  term  using  the 
relation  for  the  derivative  of  J-|(x): 


dJ1 

dx 


(2.62) 


and  obtain  the  form  identical  with  (2.13).  In  analogy  with  the 
case  of  the  straight  coastline  the  surface  boundary  condition  for 
the  temperature  is  in  analogy  with  the  case  of  the  straight  coast 
line  0  =  1  over  land  and  0  =  0  over  sea.  This  means  e(c,<t>)  =  1 
for  ?  <  r  and  0(5,0)  =  0  for  5  >  r,  where  r  is  the  scaled  radius 
of  the  island.  This  boundary  condition  leads  to  the  following  ex¬ 
pansion  coefficients  (e.g.  Charlton  1969): 

,  r/D 

0(t,O)  =  J^(Oj)  •  2  /  y  J0(o*y)dy  (2.63) 

o 

or  by  performing  the  integration 

e(l,  0)  =  2J1(oJtr/D)  rojj/D  •  { Ji  (  °  O  °  0~2  (2.64) 


The  solution  can  now  be  obtained  as  above  employing  the  same 
boundary  conditions.  Since  we  have  already  scaled  with  the  ap¬ 
propriate  intrinsic  length  scale  L,  the  flow  structure  depends 
only  on  the  scaled  radius  r.  For  large  radii  r  >>  1  the  effect 
of  curvature  must  be  negligible  and  the  results  for  all  fields 
near  the  coastline  should  be  identical  to  the  results  obtained 
for  the  case  of  a  straight  coastline.  For  radii  comparable  to 
1  or  smaller  the  maximum  radial  velocity  should  be  influenced 
by  the  effect  of  curvature  of  the  coastline.  The  result  for  the 
radial  velocity  at  5  =  r  and  n  ■  0.3  (where  the  maximum  in  u  is 
found  for  large  r)  is  shown  in  Fig.  2.18.  No  appreciable  effect 
of  the  curvature  is  found  when  r  t  1 ,  and  the  point  where  umax 
is  half  the  value  at  infinite  r  is  the  value  r  =  0.06.  Using 
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Fig.  2.18.  Radial  velocity  ( u-component)  for  the  circular 
island  at  £=radius,  and  n=0.3,  shown  for  diffent  radii. 

j  the  values  K  =  5  m2  s'"1  and  N  *  10“2  s"1  the  scale  L  becomes 

(  approximately  36  km  and  the  radius  has  to  be  less  than  ~  10  km 

|  before  the  maximum  inflow  velocity  is  influenced  by  the  finite 

size  of  the  heated  area.  The  rather  small  value  of  r  for  which 
curvature  effects  are  important  for  the  maximum  inflow  velocity 
shows  the  importance  of  small  scales.  As  discussed  in  section 
2.6  these  small  scales  are  attenuated  strongly  in  the  presence 
.j  of  even  a  small  basic  state  velocity  across  the  coast.  Also  non¬ 

linear  interactions  become  important  for  these  scales  even  with 
^  small  values  of  the  surface  temperature  perturbation.  As  a  con- 

!  sequence  we  must  expect  the  curvature  effects  to  be  of  import- 

^  ance  at  somewhat  larger  valves  of  r  in  the  real  atmosphere.  The 

analysis  could  be  extended  as  in  section  2.6  by  the  introduction 
|  of  a  basic  state  velocity  U;  however,  this  will  lead  to  a  much 

more  complicated  mathematical  problem  since  obviously  the  circu- 
!  lar  symmetry  assumed  here  will  be  violated  in  this  case,  and  the 

flow  will  become  truly  three-dimensional.  It  should  be  noted  in 
passing  that  a  rather  simple  technique  has  been  devised  for 
problems  of  this  kind  by  Scorer  (1956)  and  Olfe  and  Lee  (1971). 
i  This  latter  work  is  of  particular  interest  here  since  it  deals 

with  the  urban  heat  island  studied  by  linearized  equations  as  in 
|  the  present  investigation.  This  technique  is  based  on  the  super- 

|  i  position  of  solutions  for  two-dimensional  surface  heating  distri- 
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butions.  For  the  island  flow  discussed  here  the  appropriate  two 
dimensional  problem  is  the  flow  over  an  infinitely  long  strip  of 
land  of  finite  width  of  say  21.  This  can  be  solved  directly  as 
in  the  case  of  the  infinite  land  with  the  only  difference  that 
the  expansion  coefficients  in  the  Fourier  expansions  now  depends 
on  i  as  ( sin(  k£/M/k)  cosk?  for  the  temperature,  pressure,  and 
vertical  velocity,  and  as  -  (  sin  (  k  S./L)  sink  5  for  the  horizontal 
velocity.  The  solution  can  be  found  for  any  value  of  the  basic 
state  cross  coastal  velocity.  Noting  that  any  component  of  vel¬ 
ocity  along  the  coast  will  have  no  effect  on  the  flow  field  the 
solution  can  be  found  for  any  orientation  of  a  basic  state  vel¬ 
ocity  vector.  Superposing  solutions  for  all  orientations  of  the 
strip  of  land  keeping  constant  the  basic  state  velocity  in  an 
absolute  frame  of  reference  leads  to  a  three-dimensional  solution 
for  the  flow  over  an  isolated  island  with  the  surface  temperature 
distribution  given  by  the  superposition  as: 

<  * 

1  for  5  *  - 

9(5,0)  =  { 

2  l  l 

—  arcsin  —  for  5  >  — 
u  L5  L 

This  function  is  shown  in  Fig. 2. 19.  By  using  distributions  for 
the  two-dimensional  surface  temperature  field  other  than  the 
"top  hat"  used  here  other  three-dimensional  distributons  are 
obtained;  the  method  will  always  give  distributions  with  a  tail 
~  5-1 .  In  the  urban  heat  island  such  distributions  may  ad¬ 
equately  approximate  the  rather  smooth  surface  temperature  dis¬ 
tributions  observed,  and  in  the  study  by  Olfe  and  Lee  (1971)  the 
three-dimensional  temperature  field  over  a  city  heat  island  was 
studied  using  this  technique.  Their  model  differs  from  the  lin¬ 
earized  one  described  here  in  that  rotation  and  viscous  effects 
were  neglected  and  the  surface  forcing  was  constant  in  time 
(to  *  0).  In  this  study  we  will  not  proceed  along  these  lines 
but  note  that  the  periodically  forced  model  can  be  employed  for 
this  rather  complex  flow  problem  and  give  semi-analytical  sol¬ 
utions  for  the  perturbation  flow  with  proper  boundary  conditions. 
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Fig.  2.19.  Surface  temperature  distribution  obtained  from  the 
superposition  of  infinite  strips  of  width  2i.  (See  text). 


3.  NUMERICAL  MODEL 


3.1.  Basic  equations 


The  model  is  based  on  the  anelastic  approximation  (Ogura  and 
Philips  1962,  Gough  1969).  The  approximation  is  justified  by 
scale  analysis  of  the  terms  in  the  Navier  Stokes  equations 
treating  the  flow  variables  as  perturbations  on  a  basic  state 
of  constant  potential  temperature.  The  scale  height  of  poten¬ 
tial  temperature  is  roughly  four  times  higher  than  the  density 
scale  height,  and  the  anelastic  approximation  is  consistent  with 
the  continuity  equation  for  "deep  convection"  (Dutton  and  Fictl 
1969) s 


73  *  <°o 


(3.1) 
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where  p0  is  the  height  dependent  density  corresponding  to  the 
basic  state  of  constant  potential  temperature.  The  inelastic 
approximation  has  been  used  in  models  of  gravity  waves  and  deep 
(cumulus)  convection.  The  approximation  effectively  filters  out 
sound  wave  solutions  and  is  not  based  on  the  hydrostatic  approxi¬ 
mation.  Neuman  and  Mahrer  (1971)  argue  that  it  is  necessary  to 
include  nonhydrostatic  effects  in  their  model  of  the  sea-breeze; 
however,  as  pointed  out  by  Walsh  (1974)  their  argument  is  false, 
and  from  the  point  of  view  of  scaling  of  the  Navier  Stokes 
equations  effects  of  dynamical  pressure  is  unimportant  when  the 
aspect  ratio  of  the  flow  (typical  height/typical  width)  is 
small.  It  cannot,  however,  be  ruled  out  that  nonhydrostatic 
effects  are  important  in  parts  of  the  flow,  and  when  modelling 
in  detail  the  propagating  sea-breeze  front  it  should  be  expected 
that  nonhydrostatic  effects  come  into  play.  The  modelling  of 
such  small-scale  features  requires  a  very  high  resolution,  and 
since  most  model  formulations  require  that  model  boundaries  are 
effectively  outside  the  influence  of  the  flow,  the  number  of 
gridpoints  required  (or  equivalently  the  number  of  functions  in 
a  Fourier-type  expansion)  becomes  very  high.  A  typical  extension 
for  a  sea-breeze  flow  is  100  km  and  in  order  to  resolve  in  de¬ 
tail  a  frontal  structure  with  dimension  of  the  order  of  1  km  the 
gridsize  required  is  of  the  order  0.1  km  giving  ~  103  gridpoints 
assuming  constant  gridsize.  The  requirements  on  computer  time 
and  storage  very  easily  becomes  prohibitive.  Two  possibilities 
are  to  use  a  nonuniform  grid  (Anthes,  1978),  or  a  system  of 
nested  grids  with  different  gridsizes  (Walsh,  1974).  The  obvious 
disadvantage  is  a  tendency  to  concentrate  the  fine  structure  in 
the  region  with  a  fine  grid,  excluding  the  possibility  for  the 
front  to  propagate  away  from  the  coast  with  a  correct  speed.  Most 
published  sea-breeze  models  utilize  a  grid  with  a  gridsize  effec¬ 
tively  excluding  flow  structure  with  an  aspect  ratio  of  order 
one.  Typical  grid  size  of  a  few  kilometers  in  the  horizontal  com¬ 
pared  to  the  (possibly  correctly  resolved)  vertical  dimension  of 
the  flow  of  less  than  one  kilometer  limits  the  aspect  ratio  of 
the  resolved  flow  to  ~  0.1  or  smaller. 


p0  *  height-dependent  constant  basic  state  density  corre¬ 
sponding  to  the  basic  state  potential  temperature 


m 
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where  pQO  is  a  standard  surface  pressure  and  Hs  the  scale 
height  Cp  0  g_1 .  Here  Rg  is  the  gas  constant  and  Cp  the 
the  heat  capacity  at  constant  pressure. 


_  i 


The  it  functions  are  related  to  the  real  pressure  through  it  = 

-p  -  \vr P00 ) 


Cn  9  (P/Pnn)CP/Rg' 


counter  gradient  flux  correction  (see  below). 


Km,  Kn  are  the  turbulent  dif fusivities  for  momentum  and  heat, 
respectively,  calculated  as  discussed  in  the  next  section. 


3.2.  Turbulence  parameterization 

The  one-dimensional  version  of  the  model  results  from  assuming 
horizontal  homogeneity  as  : 


3  u  3  3  u 

—  =  f (v  —  V— )  +  —  K  — 
3 1  9  3z  m  3 z 


3  v  3  3v 

_  =  -  f(u  -  ug)  +  -  Kjn  _ 


(3.3) 


3R  3  /3R  \ 

It  "  Tz  Kn  \Tz  ~  Tc) 


The  turbulence  closure  is  obtained  by  the  utilization  of  the 
equation  for  the  turbulent  energy  in  connection  with  a  prognostic 
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with  Km  =  t/cie.  The  first  term  in  Eq.  (3.4)  is  the  shear  and 
buoyancy  production  terms  assuming  again  fluxes  to  be  pro¬ 
portional  to  the  local  gradients.  The  second  term  is  a  par¬ 
ameterization  of  the  divergence  of  the  sum  the  fluxes  <w'e'>  and 
<w'p'>/p  as  in  Shir  (1973)  and  Yamada  and  Mellor  (1975).  The 
last  term  in  Eq.  (3.4)  is  the  dissipation  term.  The  definition 
of  Km  implies  a  proportionality  between  the  mixing  length  scale 
and  the  dissipation  length  scale.  Eq.(3.4)  can  be  rewritten  as 


3e  _  es-e  3  3e 

3 1  TTcyTe  3z  ^  3  z 


with 


/cic31 


(3.6) 


(3.7) 


Neglecting  the  transport  term,  Eq.(3.6)  describes  the  relaxation 
of  turbulent  energy  towards  the  value  es,  where  production  and 
dissipation  are  in  balance,  with  a  time  scale  t  =  t/c3/e. 
Equation  (3.5)  for  the  length  scale  is  based  on  the  assumption 
that  the  length  scale  likewise  tends  to  an  equilibrium  value  ts 
with  the  same  relaxation  time  scale.  *s  is  assumed  to  be  deter¬ 
mined  from  the  boundary  layer  height  H: 

min(C4  •  H,  kz*"1),  for  z  <  H 

*s  "  {  (3-8) 

0  ,  for  z  *  H 


The  length  scale  equation  (3.5)  was  originally  proposed  by  Busch 
et.  al.  (1977).  They,  however,  combined  this  equation  with  the 
assumption  e*es,  which  under  conditions  of  weak  surface  heating, 
and  moderate  to  strong  wind  can  be  further  simplified  by  neglect¬ 
ing  the  buoyancy  production  term  in  (3.7).  Also  their  definition 
of  the  equilibrium  length  scale  was  different  and  given  as: 

l8  -  kz ♦“ 1  (l  -  ,  for  z  <  H  .  (3.9) 


Near  the  surface  both  formulations  give  is  =  kz/  in  consist¬ 
ency  with  surface  layer  similarity.  The  constants  c-|, 02,03  and 
c4  are  found  from  known  surface  relations  under  equilibrium 
conditions.  The  transport  term  in  Eq.(3.6)  can  be  neglected  for 
z/L<<  1  (e.g.  Wyngaard  and  Cote  1970),  and  we  have  e=es  and 
t  =  is;  introducing  this  into  (3.7)  and  rewriting  yields 


~  -s'  (’  -  £) 


(3.10) 


2  2 
i  .  —  1 1 


Km  5  u*  =  u*  /c1  ‘ 


^■S'O-Ee) 


(3.11 ) 


For  z/L<<  1  we  have  eg=  /c^c^u*  a  5u*  (Busch,  1973)  and  from 
Eq.  (3.11)  Cj/2C31  =  1  ,  yielding  c^O.2  and  c3=0.23/2.  The  con¬ 
stant  c4  in  Eq.  (3.8)  is  determined  by  requiring  the  length  scale 
to  scale  with  the  peak  in  the  vertical  velocity  spectrum,  which 
under  convective  conditions  is  found  to  approach  a  maximum  value 
Xm  =  1.5H  in  the  upper  part  of  the  boundary  layer  (Kaimal  et.al., 
1976),  whereas  for  z<<L  Xm  =  2z  (Kaimal  et.  al.,1972).  This  gives 
c4= 1 . 5k/2=0 . 26 .  The  constant  C2  is  set  equal  to  0.5  here;  results 
are  very  insensitive  to  the  choice  of  this  constant.  The  counter¬ 
gradient  flux  term  rc  appearing  in  Eqs.  (3.3)  and  (3.7)  is  taken 
from  Bodin  (1978)  as 


(3.12) 


-here  w**(g/T0<w ' e '>H) 1//3  is  the  convective  velocity  scale.  The 
justification  for  the  use  of  this  counter  gradient  flux  term  in 
the  convective  boundary  layer  is  discussed  in  Deardorff  (1966). 
The  height  of  the  boundary  layer  is  an  important  parameter  en¬ 
tering  into  the  dynamical  equations  aa  the  limiting  mixing  length 
under  convective  conditions, and  is  determined  from  the  tempera¬ 
ture  profile  as  the  height  of  the  intersection  of  the  temperature 
profile  with  the  line  passing  through  the  minimum  with  the  line 
passing  through  the  minimum  temperature  and  having  a  slope  equal 
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to  yc.  This  procedure  is  similar  to  that  devised  by  Bodin 
(  1978) . 

Under  neutral  and  stable  conditions  H  is  determined  as  the 
height  where  the  bulk  Richardson  number  RiB  exceeds  the  value  1 
with 


RiB 


(e  -  e0)  •  z 

(u2  +  V2) 


(3.13) 


The  value  1  for  the  "critical"  bulk  Richardson  number  is  somewhat 
higher  than  the  values  found  by  Mahrt  (1980).  It  was  chosen  by 
calibrating  against  the  observed  nocturnal  boundary  layer  height 
for  Wangara  day  33/34  (see  below).  The  flux  gradient  relation¬ 
ships  and  <f>h  are  based  on  the  expressions  from  Businger  et  al. 
(1971).  The  expressions  for  unstable  stratification  are  changed 
by  using  the  exponent  -1/3  in  both  4m  and  ^  in  contrast  to  the 
original  values  of  -1/4  and  -1/2,  and  changing  the  coefficients 
to  z/L  to  give  an  unchanged  slope  for  z/L  ♦  0,  this  leads  to  the 
expressions : 


(3.14) 


♦h 


0.74 


{ 

0.74 


+ 


1/3 


(3.15) 


Under  strongly  convective  conditions  -L  could  become  arbitrarily 
small  in  the  model,  whereas  in  the  real  atmosphere  the  convective 
motions  on  the  scale  of  the  depth  of  the  boundary  layer  will 
induce  horizontal  motions  near  the  surface,  which  in  turn  will 
generate  small-scale  turbulence.  This  mechanism  is  discussed  in 
Businger( 1973) .  The  effect  is  simulated  in  the  model  by  re¬ 
defining  the  friction  velocity  entering  into  the  definition  of  L 


as: 
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u?  =  u*  +  0.002  w* 


(3.16) 


r 


here  us  is  the  redefined  velocity  scale,  u*  the  usual  friction 
velocity  and  w*  the  convective  velocity  scale.  From  Eq.  (3.16) 
the  minimum  value  for  -L  becomes 


Sni 


in 


=  2.5 


10" 


(3.17) 


More  properly,  the  minimum  -L  should  depend  on  the  ratio  H/z0 
(Businger,  1973);  this  can  be  argued  as  suggested  by  Jensen 
(1982,  personal  communication)  as  follows:  Under  strongly  con¬ 
vective  conditions  the  geostrophic  drag  law  can  be  written  as 
(Wyngaard  et.al.,  1974) 


1 .  L) 

u  V  z„/ 


(3.18) 


with  vanishing  forcing  G  ♦  0  and  u*  *  0  the  average  wind  profile 
induced  by  the  convection  should  match  the  velocity  of  order  w* 
well  above  the  surface  layer, and  the  induced  velocity  scale  us 
is  then  determined  by  substituting  us  for  u* : 


£  ■  **  Q  M  ■  “*■  (*fc  ?) 


(3.19) 


Here  a  is  a  constant  of  order  one.  Eq.(3.19)  can  be  solved  . 
give  us=f (H/z0)w*  and  L,„in  “g(H/z0).  The  expression  gives  only 
the  scale  velocity  in  the  asymptotic  case  of  vanishing  u*  and 
Eq.  (3.16)  is  a  simple  but  somewhat  arbitrary  way  of  ensuring 
this  asymptote.  Businger( 1973)  argues  that  the  effect  can  ex¬ 
plain,  at  least  qualitatively,  why  the  asymptotic  behaviour  of 
♦h»  ♦h  ~  (“2/k)-1/^  found  from  measurements  differs  from  what  is 
expected  from  free  convection  scaling  ~  (-z/L)”1/^;  it  is 
clear,  however,  that  a  quantitative  evaluation  of  the  effect  re¬ 
quires  observation  of  H* ^  in  addition  to  surface  micrometeoro- 
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Or  measurements  of  the  low-frequency  power  spectrum  for  the 
horizontal  wind  components. 
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logical  measurements,  and  the  results  cannot  be  composited  using 
Monin-Obukhov  scaling. 


3.3.  Simulation  of  Wangara  day  33/34 

The  extensive  field  programme.  The  Wangara  Experiment  (Clarke, 
1971)  was  in  part  aimed  at  providing  good  datasets  on  the  bound¬ 
ary  layer  structure  under  rather  ideal  conditions,  and  the  data 


Fig.  3.1.  Wind  speed  profiles  at  14  h  day  33  and  02  h  day  34 

of  the  Wangara  experiment. ( - ):  Model  results, ( *- *-) , 

(o-o-o):  Observations.  (A):  Geostrophic  wind  speed  at  14  h, 
(A):  Geostrophic  wind  speed  at  02  h. 


65 


obtained  naturally  lends  themselves  to  the  evaluation  of  the  qual¬ 
ity  of  model  simulations.  For  a  number  of  boundary  layer  models 
day  33/34  of  the  experiment  has  been  used  for  this  purpose.  The 
main  reason  .for  choosing  these  two  days  is  the  relative  absence 
of  large  synoptic  scale  or  mesoscale  phenomena,  which  would 
seriously  complicate  the  interpretation  of  the  boundary  layer 
processes.  For  the  present  purpose  of  the  design  of  a  mesoscale 
model  particularly  aimed  at  simulation  of  the  sea-breeze,  the 
diurnal  variation  of  the  boundary  layer  structure  in  response  to 
the  surface  heating,  and  in  particular  the  response  of  the  wind 
profile  is  of  interest.  The  wind  profile  was  measured  at  the 
experiment  using  theodolite  tracking  of  balloons,  which  were 
filled  to  rise  at  a  constant,  nominal  rate  enabling  the  wind 
profile  to  be  determined  up  to  2  km.  The  initial  data  for  the 
model  was  the  observed  wind  and  temperature  profile  at  06  hr  day 
33.  The  surface  heating  function  was  prescribed  as  pure  sinus¬ 
oidal  during  the  day  with  maximum  heating  <w'  e  '>max=0 .  1 5  Kras-1 
and  following  the  observed  net  radiation.  The  surface  radiational 
cooling  at  night  was  simulated  by  prescribing  an  exponential  de¬ 
cay  of  the  surface  temperature,  with  decay  rate  and  asymptotic 
value  determined  from  the  observed  surface  temperature  at  night 
between  days  33  and  34.  This  ensures  that  the  night  temperature 
at  the  surface  in  the  model  closely  follows  the  observed  tempepe- 
ture  provided  that  the  surface  temperature  at  sunset  is  correct. 
The  geostrophic  wind  was  specified  as  varying  linearly  in  time 
and  varying  linearly  between  1  km  and  2  km,  using  the  3  hourly 
surface  geostrophic  winds  from  the  synoptic  network  given  in  the 
dataset* ) . 

The  model  was  found  to  predict  accurately  the  development  of  the 
daytime  mixed  layer  height  and  temperature  profile,  and  also  the 


*)  The  surface  geostrophic  wind  was  found  by  fitting  a  straight 
line  to  the  observed  values  of  each  of  the  two  components  in  the 
interval  between  03  day  33  and  03  day  34 .The  thermal  wind  was 
observed  only  at  12  hourly  intervals  and  the  values  used  here 
were  simply  specified  by  the  line  connecting  the  values  at  09 
day  33  and  21  day  33. 
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development  of  the  night  mixed  layer  height.  Figure  3.1  shows 
the  predicted  and  observed  windspeed  profiles  of  day  33  hour  14, 
where  the  mixed  layer  is  well  developed,  and  at  hour  02  of  day 
34,  where  a  nocturnal  jet  has  developed.  The  model  shows  good 
agreement  with  the  observed  wind  profiles  in  particular  in  view 
of  the  uncertainty  with  which  the  geostrophic  wind  is  measured. 
The  most  pronounced  difference  is  the  observed  decrease  of 
windspeed  with  height  above  1200  m,  a  feature  which  is  consistent 
through  the  night.  The  thermal  wind  specified  in  the  model  was, 
as  mentioned,  computed  from  the  synoptic  radiosonde  network, 
using  six  stations  with  separations  of  around  500  km,  and  with 
three  of  them  located  at  the  coast  (see  Clarke  (1971)  for  de¬ 
tails).  The  observed  decrease  of  windspeed  above  1200  m  could 
therefore  possibly  be  explained  by  baroclinicity  not  resolved  by 
the  radiosonde  network.  A  particular  feature  prominent  in  the 
data  is  the  development  of  a  strong  nocturnal  jet,  this  is  well 
reproduced  in  the  model,  with  regard  to  the  level  of  maximum 
windspeed,  time  of  maximum,  and  magnitude. 

3.4.  Two  dimensional  model 

The  set  of  dynamical  equations  (3.2)  for  the  mean  flow,  and  the 
turbulence  parameterization  described  above  form  the  basis  of  the 
model.  Added  in  Eqs.  (3.4)  and  (3.5)  are  advection  by  the  mean 
flow  of  the  turbulent  energy  and  the  length  scale,  respectively. 
The  numerical  solution  of  the  system  of  equations  is  based  on 
finite  difference  methods.  To  some  extent  the  different  physical 
mechanisms  are  treated  in  sequence  using  time  splitting  (Nesinger 
and  Arakawa, 1 976 ) .  The  first  step  consists  of  the  calculation  of 
the  hydrostatic  pressure  by  integration  from  the  model  top  of 
the  hydrostatic  equation,  assuming  a  linear  temperature  profile 
between  gridlevels.  The  advection  terms  are  calculated  using 
centered  differences  in  time  and  space,  and  the  result  after  the 
first  step  are  fields  based  on  the  inclusion  of  advection,  hydro¬ 
static  pressure  gradient,  and  coriolis  terms  only.  The  next  step 
consists  in  the  inclusion  of  the  diffusion  terms  by  first  comput¬ 
ing  the  forcing  by  surface  fluxes.  The  numerical  method  used  for 
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the  diffusion  terms  is  the  finite-element  method  with  implicit 
time  differencing.  The  diffusion  equation  can  be  written: 


3f  _  3  3 f 

3t  Iz  Iz 


(3.20) 


assuming  both  f(z)  and  K(z)  are  expanded  in  a  series  of  some 
basic  functions  4>j_(z),viz.  f(z)=f^  ^(z)*)  and  K(z)=Ki  ^(z). 

Then  a  set  of  equations  can  be  obtained  for  the  temporal  develop¬ 
ment  of  the  expansion  coefficient  fk  by  inserting  the  expansions 
into  Eq. ( 3 . 20 ) ,  multiplying  by  ^  and  integrating  over  z  to  yield 

3fi 

Aik  —  =  K*  fi  BUk  <3*21> 

where  the  matrices  A  and  B  are  constant  and  given  by 

Aik  =  /  ♦  i  ♦k  dz 
z 

(3.22) 

Buk  s;*k^4*idz 

These  expressions  are  quite  general  for  the  spectral  method  for 
the  solution  of  the  diffusion  equation,  and  Eq.(3.20)  can  easily 
be  shown  to  give  the  least-square  deviation  for  the  expansion 
coefficients  provided  the  expansions  of  f  and  K  are  correct. 

With  orthonormal  basic  functions,  A  would  become  the  identity 
matrix  and  no  matrix  inversion  would  be  neccesary,  but  B  would 
generally  be  complicated.  In  the  finite  element  method  the  basic 
functions  are  chosen  so  as  to  make  both  A  and  B  sparse.  In  the 
present  model  the  $'s  are  chosen  to  be  the  "chapeau"  functions. 
The  k'th  chapeau  function  is  zero  everywhere  except  in  the  inter¬ 
val  form  the  (k-l)'th  to  the  (k+1)'th  gridpoint,  it  is  piecewise 
linear  ,and  takes  the  value  of  1  at  the  k'th  gridpoint.  With  a 
function  given  as  gridpoint  values  the  expansion  coefficients 
for  the  function  on  this  grid  simply  become  the  gridpoint  values 
themselves,  and  the  expansion  is  correct  only  if  the  function  is 


*)  Here  summation  over  repeated  indices  is  assumed. 
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piecewise  linear.  With  this  choice  for  the  matrix  A  becomes 
tridiagonal,  with  main  diagonal 

Akk  =  V3  ( izk-1+Azk>  and  Ak-1  ,k  =  1/6Azk_i ,Ak+i fk  =  1/6azk. 

The  elements  of  B  are  given  by  Bk-1 ,k-1 ,k  =  Bk,k-1 ,k  = 

"  Bk-1,k,k  =  1  / ( 2Az)C_i  ) ,  Bk+1/k+1fk  =  Bkfk+1fk 
=  ~Bk+1 ,k,k  =  1/(2Azk),  and  Bk>kfk  =  Bk_i fkfk  +  Bk+i ,k,k  • 

All  other  terms  are  zero.  The  implicit  scheme  for  the  time 
differencing  is  a  special  case  of  the  Cranck-Nicholson  scheme, 
employed  here  for  the  extrapolation  from  timestep  n  to  n+1  : 

Aik  (f[n+1)  “  fin))  =  Kin)  fin+1)  Biik  At  (3.23) 

The  next  step  consists  in  including  the  effect  of  the  dynamical 
pressure  term.  After  this  step  the  fields  at  timestep  n+1  are 
given  by  the  result  after  the  diffusion  step  plus  the  tendency 
from  the  x-terms: 


n+1  =  w*  -  (?s  At) 


where  u*,w*  are  the  fields  obtained  after  the  diffusion  step. 
Applying  the  divergence  operator  given  as  D  *  (3x,po_l3z  po^ 
to  the  (u,w)  in  Eq. (3.24)  results  in  the  equation 


3  3  n  „n+1  *  3  U*  +  p""^  3  P  W* 

X  *0  °ZW0W  X  v0  °zy0w 


-3^«-p^3o3ir 
XX  Mo  zMo  z 


(3.25) 


where  ir  has  been  redefined  by  letting  it  absorb  the  at.  The 
continuity  equation  requires  that  the  left-hand  side  vanishes, 
and  we  obtain  a  differential  equation  for  w.  The  finite  differ¬ 
ence  operators  corresponding  to  3X  and  3Z  are  centered  dif¬ 
ferences.  With  D*  *  D(u*,w*)  calculated  *  is  obtained  from  the 
the  solution  of  this  Poisson-type  differential  equation.  The 
boundary  conditions  for  *  are  3Z  *  at  the  lowest  level  and  at 
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the  top,  and  *  =  0  at  the  lateral  boundaries.  The  solution  is 
obtained  by  a  direct  method,  which  consists  in  first  performing 
a  horizontal  Fourier  transform  at  each  level.  The  system  re¬ 
sulting  for  each  Fourier  mode  can  be  written  as 

— 1— -  ( cos ( 2 Axk )- 1 )  +  P~J  3zp03z l'k  =  Dk  (3.26) 

l2Ax2  j 

The  first  term  in  the  parenthesis  is  the  response  function  for 
the  horizontal  difference  operator  3XX,  and  d£  is  the  k'th 
Fourier  mode  of  D*.  The  boundary  conditions  for  u  and  w  (see 
below)  gives  D*  zero  at  the  lateral  boundaries  and  the  Fourier 
transformation  can  be  written  involving  only  sine  terms.  With  a 
centred  difference  operator  for  3Z  Eq.(3.26)  results  in  a  sys¬ 
tem  of  linear  equations  with  a  pentadiagonal  coefficient  matrix, 
which  can  be  solved  by  standard  methods.  The  final  step  then 
consists  in  applying  the  inverse  sine  transform  to  the  expansion 
rrk  at  each  level.  To  improve  the  speed  in  the  computations  the 
operator  p0-13z  p03z  e<3-(3.26)  was  written  as  (p0-13zp0)3z 
+  3  zz>  and  with  p0  -13zp0  computed  exactly  from  the  analyti¬ 
cal  expression  for  p0  at  each  level.  Using  a  finite  element 
method  similar  to  the  method  described  above,  the  resultant  lin¬ 
ear  equations  have  a  tridiagonal  coefficient  matrix  which  is  con¬ 
stant.  This  system  is  conveniently  solved  by  using  the  the 
method  described  by  Ahlberg  (1977)  performing  the  first  decom¬ 
position  step  once  and  for  all. 


3.4.1.  Boundary  conditions 

The  computational  domain  is  bounded  at  the  top  by  a  rigid  lid 
where  the  initial  values  for  u,v  and  8  are  assumed  to  be  main¬ 
tained  and  where  "h  =  '•z*  “  w  *  0.  At  the  surface  the  usual 
no-slip  conditions  are  assumed,  i.e.  u  =  v  *  w  =  0,  and  3Z*  =  0. 
The  surface  boundary  conditions  for  the  temperature  are  that  over 
water  the  initial  surface  temperature  is  maintained,  whereas  over 
land  the  boundary  conditions  are 


during  the  day:  <w'8'>0  ■  <w'0,>max  8*n  (j“ 
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at  night:  T0  -  +  (Tgunset--Tinitial )  exp(-t/tc) 

similar  to  the  boundary  conditions  used  in  the  simulation  of 
Wangara  day  33/34.  These  boundary  conditions  reflect  the  physi¬ 
cal  mechanisms  of  a  daytime  heated  boundary  layer  with  heating 
governed  by  the  solar  input  at  the  surface,  and  a  night-time 
boundary  layer,  where  the  ground  surface  temperature  is  primarily 
governed  by  direct  long  wave  radiation  to  space.  The  lateral 
boundary  conditions  are  3xu  =  3xv  =  3X0  =  w  =  u  =  0.  The  first 
ten  meters  above  the  surface  are  assumed  to  obey  surface  layer 
similarity  theory,  and  the  actual  boundary  conditions  entering 
into  the  prognostic  equations  are  the  fluxes  of  momentum  and  heat 
at  the  first  level  at  10  m.  This  procedure  is  commonly  used  in 
boundary  layer  models  as  in  the  model  by  Busch  et.  al.(1977).  The 
surface  boundary  condition  for  w  is  used  as  w(10m)  =  0,  thus  as¬ 
suming  negligible  convergence  in  the  10m  layer. 


3.4.2.  Horizontal  diffusion  and  filtering 

The  model  equations  (3.2)  include  no  horizontal  diffusion  terms. 
As  discussed  in  section  2.1  the  effect  on  the  mean  flow  struc¬ 
ture  from  horizontal  diffusion  should  be  unimportant  if  the  flow 
aspect  ratio  is  much  smaller  than  unity.  In  a  numerical  model  the 
combined  effect  of  truncation  errors  in  the  finite  difference 
operators,  and  aliasing  introduced  by  the  nonlinear  terms,  will 
result  in  large  errors  if  no  precaution  is  taken  to  damp  the 
small-scale  features.  In  the  present  model,  as  in  most  mesoscale 
models  reported  in  the  literature,  the  resolution  is  insufficient 
to  accurately  resolve  all  the  physical  features  which  we  want  to 
model.  We  must  therefore  rely  on  the  asumption  that  the  smallest 
resolvable  scales  have  only  a  negligible  influence  on  the  behav¬ 
iour  of  the  larger  scales  of  motion.  The  relative  success  of 
large-scale  weather  prediction  models,  in  spite  of  similar  res¬ 
olution  problems  with,  for  example,  the  proper  resolution  of 
fronts,  suggests  that  this  assumption  is  often  justified,  and 
the  results  with  the  present  model  in  the  simulation  of  the 
propagation  of  the  seabreeze  front  described  below  also  suggests 
that  this  is  a  minor  problem.  The  one-dimensional  advection-dif- 
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fusion  equation  with  a  constant  advection  velocity  U  and  dif- 
fusivity  K  can  be  written 


ax 

at 


a  y 

u  _ _  +  K 

3x 


(3.27) 


With  a  finite  difference  approximation  of  this  equation,  trunc¬ 
ation  errors  in  the  advection  term  will  produce  small-scale 
noise  at  a  rate  which  must  be  proportional  to  U/Ax.  Similarly, 
the  destruction  of  small  scales  by  the  diffusion  term  must  be 
proportional  to  K/Ax2  .  If  we  define  a  grid  Reynolds  number 
as  the  ratio  of  these  two  numbers: 


U  Ax 

Re  =  _ 

K 


(3.28) 


we  can  control  the  amplitude  of  the  small  scales  by  specifying  a 
maximum  value  of  Re,  that  is  letting  K  depend  on  u  through  (3.28) 
with  Re  fixed.  The  linear  advection-dif fusion  equation  (3.27) 
will  not  give  rise  to  aliasing,  but  in  the  more  general  form  in 
the  model  equations  (3.2)  aliasing  will  be  present  in  the  finite 
treatment  of,  for  example,  the  u3u/ax-term,  and  also  in  the  other 
nonlinear  terms.  In  the  present  model  the  neccesary  horizontal 
smoothing  is  accomplished  by  specifying  Re  as  constant  throughout 
the  integration,  and  using  the  maximum  value  of  u(x,z,t)  at  every 
timestep  as  U  in  (3.28),  giving  the  value  of  K  used  as  the  hori¬ 
zontal  diffusivity.  The  diffusion  term  is  computed,  using  a  3- 
point  formula  for  3XX  and  an  explicit  forward  time  extrapolation, 
as  the  final  step  in  the  sequence  described  in  section  (3.4). 


The  above  consideration  s  also  applies  to  finite  difference  oper¬ 
ators  acting  in  the  vertical  direction.  In  the  model  a  minimum 
diffusivity  defined  for  each  vertical  column  of  grid  points  is 
adopted,  using  (3.28)  with  UAx  replaced  by  &Zk)max.  The 

turbulent  dif fusivities  from  the  parameterization  described  in 
Section  (3.2)  are  therefore  limited  from  below  by  this  minimum 
value  before  being  applied  in  the  model.  To  avoid  the  accumu¬ 
lation  of  noise  near  the  lateral  boundaries,  an  increased  hori¬ 
zontal  diffusion  is  applied  in  a  zone  near  each  of  the  lateral 
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boundaries.  This  is  accomplished  by  the  application  of  a  3-point 
filter  (Shapiro,  1970)  at  every  timestep  in  these  diffusive  zones 
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3.5.  Comparision  with  analytical  model 

The  results  from  the  linearized  model  described  in  the  previous 
chapter  should  be  approximately  reproduced  by  the  numerical  model 
if  the  surface  forcing  is  so  weak  that  the  nonlinear  advection 
terms  are  negligible,  and  the  turbulent  dif fusivities  are  exter¬ 
nally  specified  as  in  the  linear  model.  Comparisons  between  the 
models  then  yield  an  independent  check  on  the  numerical  pro¬ 
cedures  with  regard  to  accuracy,  stability,  and  code  errors.  The 
stability  and  accuracy  of  the  individual  finite-difference  oper¬ 
ators  in  the  numerical  model  can  be  performed  by  standard 
methods,  and  this  kind  of  analysis  is  important  in  designing  the 
model,  but  the  adequateness  of  the  model  as  a  whole  can  probably 
be  asserted  only  by  using  known  solutions.  The  layered  linear 
model  enables  a  realistic  vertical  profile  of  the  diffusivity 
and  stratification  in  the  boundary  layer.  This  is  important  for 
the  detailed  vertical  structure  of  the  mean  profiles,  and  specifi 
cally  for  the  maximum  wind  velocities  as  discussed  in  Section 
2.8.  Furthermore,  the  most  important  physical  mechanism  in  the 
sea-breeze  flow  is  the  vertical  diffusive  transport  of  heat  and 
momentum;  therefore,  it  is  important  that  transport  in  particu¬ 
lar  be  tested  as  realistically  as  possible.  Figure  3.2a  shows 
examples  of  such  comparisons  in  the  case  where  the  dif fusivities 
(K^*Km*K)  are  assumed  to  be  closely  parabolic  up  to  the  top  of 
the  boundary  layer,  and  constant  above  the  layer  as  shown  in 
Fig.  3.2a.  In  addition  to  the  diffusivity  profile,  this  figure 
shows  the  computed  analytical  u-profile  at  the  time  of  maximum 
surface  temperature  difference  together  with  the  corresponding 
result  from  the  numerical  model.  Here  the  surface  forcing  is 
assumed  harmonic  also  in  the  horizontal  corresponding  to  the 
inclusion  of  only  one  wavenumber  in  the  analytical  model  (Sec¬ 
tion  2.2);  the  wavelength  is  7  gridunits  in  the  numerical  model 
corresponding  to  28  km,  and  the  Brunt-Vaisala  frequency  is 
N-1/150  s'1 .  Figures  3.bd  show  the  u,v,  and  w  components,  using 
the  same  parameters  except  that  the  stratification  below  the 
816  m  level  is  reduced  by  a  factor  of  five. 


Pig.  3.2b.  As  fig.  3.2a,  but  in  the  case  with  a  jump 
stratification  at  816  m  with  (^“ISO  s  and  N^*0.2*N2 
(See  text). 


The  profiles  are  shown  in  ms-1  per  deg  in  the  surface  temperature 
amplitude.  Because  of  the  vertical  variation  of  the  parameters, 
no  simple  scaling  exists.  Using  the  stratification  below  816  m 
and  the  maximum  diffusivity  gives  a  length  scale  of  30  km  in  the 
linear  model,  and  thus  the  nondimensional  wavenumber  is  close  to 
unity.  This  means  that  all  terms  in  the  linearized  equations  are 
of  comparable  importance.  The  u  and  w  profiles  show  very  good 
agreement  between  the  numerical  and  analytical  solutions.  The 
numerical  model  initially  assumes  zero  wind  conditions  at  the 
start  of  the  heating  period,  whereas  the  analytical  model  assumes 
periodic  solutions,  i.e.  any  transient  solutions  are  neglected, 
and  consequently  differences  should  be  expected  due  to  the 
appearance  of  inertial-diffusive  oscillations.  In  this  case  the 
numerical  model  was  integrated  for  one  period,  and  during  this 
time  the  differences  between  the  two  models  show  a  periodic 
variation  with  nearly  constant  amplitude.  The  contributions  from 
these  oscillations  are  most  pronounced  in  the  v-profile  as  should 
be  expected,  since  this  component  is  forced  only  through  the 
coriolis  term,  whereas  the  u-component  is  forced  more  directly 
through  the  pressure  gradient  term.  In  all  profiles  the  differ¬ 
ences  are  largest  at  the  higher  levels  in  the  boundary  layer, 
partly  reflecting  the  lowered  resolution  in  the  numerical  model 
there.  The  comparisions  show  that  the  the  numerical  schemes  are 
well  behaved;  in  particular  that  the  time  splitting,  and  the 
sequential  "corrector"  method  used  for  the  calculation  of  the 
dynamical  pressure  term  are  not  introducing  systematic  or  accumu¬ 
lating  numerical  errors. 

3.6.  Sea-breeze  simulations 

The  model  has  been  integrated  with  a  number  of  different  speci¬ 
fications  of  the  parameters  pertaining  to  the  synoptic  or  large- 
scale  conditions.  The  aim  has  been  to  compare  the  model  results 
with  the  qualitative  response  of  the  sea-breeze  flow  to  vari¬ 
ations  in  these  parameters.  A  large  amount  of  observational 
studies  on  the  behaviour  of  the  land-  and  sea-breezes  has  been 
reported  in  the  literature  (see,  e.g.  Schroeder  et  al,  1962  for 
references),  and  thus  some  general  statements  can  be  made  re- 
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garding  the  interdependence  of  the  synoptic  situation  and  the 
sea-breeze  flow  (as  summarized  in  Schroeder  et  al,  1962).  Because 
of  the  number  of  parameters,  which  can  be  varied  in  the  model, 
and  the  need  to  maintain  physical  realism,  a  complete  study  of 
model  response  is  prohibitive.  As  a  consequence  variation  of 
parameters  was  performed  by  using  a  reference  set  of  parameters, 
and  varying  only  one  parameter  at  a  time,  keeping  all  other  at 
their  reference  value.  The  parameters  chosen  is  given  by  the  fol¬ 
lowing  set  with  the  reference  values  listed  first,  and  the  other 
values  chosen  are  listed  in  parentheses: 


Parameters  varied: 


Maximum  surface  heatflux  over  land:  200  (100,400)  Wnf 
Basic  state  surface  geostrophic  wind 

ug  :  0  (  -3,  3  )ms“ 
v_  :  0  (  -3 ,  3  )ms” 


Basic  state  thermal  wind 


x-component  Sx  ,  below  1500  m  :  0  (  -2 ,  2  )ms-z/km 

above  1500  m  :  0 

y- component  Sy  :  0 

Basic  state  stratification  N  :0.01  (0 .02 ,0 .005 )s-1 

Grid  Reynolds  number  Re  :  2  (  5  ) 

Decay  time  for  surface  temperature  excess  over  land  at  night 
(e-folding  time)  t^  :  4.2  (16.7) 


Parameters  kept  constant: 


Heating  period 
Coriolis  parameter 
Sea  surface  temperature 
Initial  surface  temperature 
Surface  roughness  length 
over  land 
over  sea 


1/2  tjj:  12  hours 


f  :  1 .21 *10 


:  283 


:  0.01  mm 


-4  „-1 


I 


79 


Horizontal  grid 

128  points  with  a  seperation  of  2  km,  coast  at  gridpoint 
no  65.  x=0  at  coast  positive  inland. 

Diffusive  zone  (sect.  3.4.2)  at  outermost  16  gridpoints 
near  lateral  boundaries. 

Vertical  grid 

21  levels  with  increasing  separation 

z-|  =  10m  ,  Az i  =  30m  ,  Az^+1  =  Az^  *1 . 23  (  Az^  =  z ^+1 -z ^ ) 

resulting  in  Ztop=z21=  8074  m. 

Timestep 

0.8  times  max  allowable  given  as  N-1 ,  and  less  than  2  min. 

Table  3.1  shows  some  of  the  model  results  after  integrating  the 
model  for  8  hours.  The  initial  data  was  in  every  case  a  tempera¬ 
ture  field  with  constant  Brunt-Vaisala  frequency  N,  and  a  wind 
field  given  by  the  geostrophic  values  at  all  levels.  At  8  hours, 
which  is  two  hours  after  maximum  heating  over  land,  the  sea- 
breeze  is  well  developed,  and  the  maximum  of  the  onshore  wind 
component  u  at  the  coast  is  found  near  this  time  in  the  case  of 
no  synoptic  forcing.  In  the  table  each  column  gives  results  from 
choosing  reference  values  of  the  parameters  except  for  the  par¬ 
ameter  given  at  the  top.  The  first  row  shows  the  u-component  at 
the  height  of  the  maximum,  which  in  all  cases  is  near  100  m  with 
a  very  small  vertical  gradient  below  this  level.  Similarly,  the 
value  of  max(-v)  at  the  coast  is  shown  in  the  second  row.  The 
next  four  rows  show  the  maximum  u-  and  v-components,  and  the 
distance  from  the  coast  where  they  occur.  In  the  case  of  no 
large-scale  forcing  (first  five  columns),  the  maximum  winds  are 
found  near  the  coast,  and  the  maximum  values  are  close  to  those 
at  the  coast.  With  onshore  geostrophic  wind  (Ug*3ms_1),  the  maxi¬ 
mum  wind  is  displaced  inland,  and  the  strength  of  the  onshore 
component  relative  to  the  background  value  (u-Ug)  is  reduced 
roughly  by  a  factor  of  two  compared  with  the  reference  value. 

In  contrast  we  see  that  with  the  geostrophic  wind  directed  off¬ 
shore  (Ug*-3ms-1),  the  value  of  umaJ{-Ug  is  increased  by  approxi¬ 
mately  50%.  This  strong  influence  on  the  flow  from  even  weak 
synoptic  forcing  is  also  evident  from  the  differences  in  maximum 
vertical  velocity.  In  the  case  of  the  synoptic  forcing  acting  in 
the  same  direction  as  the  sea-breeze,  the  maximum  vertical  vel- 
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vertical  velocity  is  only  one-third  of  the  reference  value.  With 
Ug=-3ms_1  the  maximum  w  is  nealyr  a  factor  of  three  larger  than 
the  reference  value.  These  results  are  in  accordance  with  obser¬ 
vations  (Schroeder  et  al,  1962  ).  The  row  denoted  AT^O)  gives  the 
surface  temperature  rise  relative  to  the  initial  value  ( =  sea 
surface  temperature)  at  the  gridpoint  farthest  inland,  where  con¬ 
ditions  are  nearly  homogeneous  in  the  horizontal.  The  boundary 
layer  height  at  the  same  point  is  shown  in  the  last  row.  The  sur¬ 
face  temperature  is  only  slightly  reduced  in  the  presence  of  a 
large  scale  wind,  and  this  decrease  reflects  the  decrease  of  the 
superadiabatic  temperature  gradient  near  the  surface,  which  can 
be  seen  by  comparing  it  with  AToo(IOm),  which  is  the  temperature 
rise  at  the  same  point  at  the  10-m  level.  The  frontal  position 
is  defined  as  the  point  with  the  steepest  gradient  of  u  at  the 
10-m  level. 

j  The  position  of  the  front,  which  at  this  hour  is  weak  except  when 

|  there  is  adverse  large  scale-flow  present,  is  slightly  behind  the 

point  with  maximum  vertical  velocity,  as  expected.  The  greater 
part  of  the  horizontal  convergence  is  confined  to  the  lowest  few 
hundred  meters  above  ground.  Figure  (3.2)  shows  the  variation  of 
the  u-profile,  at  the  point  of  maximum  value,  for  the  different 
s  values  of  surface  heatflux  over  land.  The  maximum  values  are 

listed  in  Table  3.1,  and  we  see  that  increasing  the  heatflux  by 
a  factor  of  four  from  100  WnT2  to  400  Wm~2  leads  roughly  to  a 
doubling  of  both  the  maximum  windspeed  and  depth  of  the  inflow 
layer.  We  can  compare  these  with  the  results  obtained  in  section 
2.6  (Eq.  2.53),  which  can  be  interpreted  as  the  relation  that  the 
maximum  speed  umax  is  proportional  to  the  square  root  of  the  sur- 
I  face  heatflux.  This  relation  is  followed  by  the  model  to  within 

10%.  This  square  root  dependence  was  also  found  by  Pearson  (1973). 
In  contrast  to  this  result  we  see  a  substantial  dependence  of 
umax  on  the  stratification  N.  In  the  linear  model  umax  is  pro¬ 
portional  to  N-1,  whereas  Eq.  (2.53)  predicts  umax  independent 
of  N.  The  numerical  model  thus  gives  a  value  in  between  these 
two  extremes.  Table  3.2  lists  model  results  after  12  hours  of 
integration,  at  the  end  of  the  heating  period.  At  the  coastline 
the  u-component  has  decreased  compared  with  the  values  at  8 
hours,  but  the  v-component  has  simultaneously  increased  as  a  re- 
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suit  of  Coriolis  turning,  and  the  scalar  windspeed  has  increased. 
There  is  one  notable  exception,  namely  the  case  with  a  strong 
initial  stratification  N=0.02  s-1 .  In  this  case  the  v-component 
is  relatively  reduced  both  at  8  hours  and  at  12  hours.  This  re¬ 
duction  of  the  Coriolis  turning,  or  increase  of  the  cross-isobar 
angle  of  the  flow,  is  consistent  with  the  results  obtained  in 


(ms"') 

Fig.  3.3.  u-profiles  at  t»12  h  from  the  points  with  maximum 
values  near  the  ground.  Unit  on  vertical  axis  is  1  km. 

( - ).  Ref.  with  Hq*200  Wm“2  (profile  at  x-36  km), 

(-*-•-):  Ho«400  Wm“2  (50  km),  ( - ):  H0«  100  Wm-2  (22  km). 
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section  2.8.  The  point  with  maximum  u  and  v  has  moved  inland  and 
the  two  maxima  have  risen.  The  surface  front  has  likewise  moved 
further  inland,  and  has  become  stronger  as  is  evident  from  the 
large  increase  in  vertical  velocity.  The  last  row  in  Table  3.2 
shows  the  ratio  of  the  average  maximum  u-component  to  the  average 
speed  of  the  front  between  hours  8  and  12.  We  see  that  this  value 
is  close  to  2  in  all  cases  where  there  is  no  large-scale  flow  in 
the  x-direction.  This  value  agrees  well  with  the  value  of  ~  2  re¬ 
ported  by  Simpson  et  al  (1977),  and  the  value  2.2  found  by  Kimura 
and  Eguchi  (1978).  The  strong  dependence  on  the  large-scale  flow, 
which  is  evident  from  the  table,  suggests  that  this  agreement  may 
be  fortuituous,  however.  The  values  obtained  from  the  observ¬ 
ational  studies  are  based  on  the  average  of  a  large  number  of 
observations,  and  we  could  expect  a  dependency  from  the  time  of 
day  as  also  suggested  by  Physick  (1980).  Table  3.3  gives  the  rate 
of  advence  of  the  front  relative  to  the  maximum  u-component  for 
different  hours  in  the  reference  case,  and  for  the  case  with 
larger  heatflux. 


Table  3.3.  Progression  of  the  surface  sea-breeze  front. 

A:  reference  case  .  B:  Heatflux  »  400  Wm'2, 

In  the  last  case  the  model  was  run  with  the  double  number 
of  gridpoints  in  the  x-direction. 


A 


hour 

5 

7 

9 

1 1 

13 

15 

17 
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7 

19  29  43  63 

79 

91 

Hmax 

ms-1  s 

2.5 
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2.2 

1.8 
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2.3 

1.8 

1.3 
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B 

hour 

5 
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17 
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17 

30  43  61  89 
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5.3 

4.6 

4.0 

ms-1  s 

1.7 

1.8 

1.8 

2.5 

2.5 

3.1 

2.8 

“max/Ufront  * 

2.2 

2.4 

2.5 

2.0 

2. 1 

1.5 

1.4 

We  see  that  in  the  reference  case  the  ratio  umax/ufront  is  rather 
variable  with  a  value  of  ~  2  only  between  hour  8  and  12  corre¬ 
sponding  to  the  afternoon.  After  model  sunset  the  front  accele¬ 
rates  resulting  in  a  decrease  of  the  ratio.  The  value  at  5  hours 
is  not  very  well  defined  because  the  front  is  very  weak  at  this 
hour.  Figure  3.4  shows  the  u-profiles  at  12  hours  for  the  three 
values  of  stratification  at  the  points  where  u  is  largest  near 
the  surface.  We  see  that  the  profiles  are  nearly  identical  if  a 
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Fig.  3.4.  u-profiles  at  t*  12  h  and  x*32  km.  ( - )s  Ref 

with  N-0.01  a-1,  ( - );  N-0.02  s_1 ,  N-0.005  s' 


ratios  of  /K  of  1:0. 5: 1.5.  These  ratios  are  close  to  the  values 
which  we  would  estimate  from  the  figure.  Figure  3.5  shows  the 
profiles  of  u-ug  selected  as  for  Fig.  3.4  but  for  the  differ¬ 
ent  values  of  the  large-scale  wind  field.  The  u-profiles  relative 
to  the  large-scale  values  are  seen  to  be  relatively  insensitive 
to  the  value  of  the  latter,  except  in  the  case  where  they  add  to 
each  other,  in  which  case  a  large  reduction  of  umax  results.  It 
should  be  noted  that  the  profiles  are  not  taken  at  the  same  x- 
coordinate  but  at  the  points  listed  in  the  legend.  Figure  3.6 
shows  the  u,v,  and  w  fields  at  8  hours  in  more  details.  The  maxi¬ 
mum  inflow  is  concentrated  near  the  surface  as  is  also  evident 
in  the  profiles  in  Figs.  3. 3-3. 5  (at  t=12  h).  The  asymmetry  in 
both  u  and  v  is  small  only  at  this  time,  but  in  w  there  is  a 
more  pronounced  difference  between  the  up-  and  downdraught  re¬ 
gions.  The  convective  boundary  layer  over  land  slopes  upward  fol¬ 
lowing  the  zero  line  in  the  w-plot  until  the  maximum  height  of 
approximately  1.6  km  is  reached;  from  here  it  has  a  nearly  con¬ 
stant  depth  as  one  moves  inland.  The  downdraught  is  thus  con¬ 
centrated  in  a  rather  thin  layer  on  the  top  of  the  boundary 
layer  with  a  wide  extension  seaward  with  only  weak  subsidence. 

The  same  velocity  fields  6  hours  later  (2  hours  after  model  sun¬ 
set)  are  shown  in  Fig.  3.7.  At  this  time  a  stronger  front  has 
been  established,  and  has  moved  inland  as  already  described 
above.  The  definition  of  the  frontal  position  used  above,  as  the 
point  with  largest  gradient  in  u  at  10  m,  is  seen  to  be  uncriti¬ 
cal,  since  the  isolines  in  both  u  and  v  are  nearly  vertical  near 
the  surface  and  essentially  coinciding.  At  the  coastline  the 
surface  wind  is  nearly  parallel  to  the  coast  at  this  time  with 
a  comparatively  large  counter  flow  between  1  and  2  km.  The  maxi¬ 
mum  value  of  u  has  decreased  relative  to  that  at  12  hours,  but 
we  see  on  the  other  hand  that  the  maximum  value  of  w  has  actu¬ 
ally  increased  indicating  a  strengthening  of  the  front.  Compared 
to  the  situation  at  8  h,  the  w  field  at  14  h  shows  a  downdraught 
with  a  strength  much  smaller  than  the  updraught,  but  again  lo¬ 
cated  above  and  slightly  seaward  of  the  updraught  zone.  At  this 
hour  the  land  surface  has  cooled  and  a  shallow  stable  surface 
layer  has  formed.  The  shear  zone  behind  the  surface  front  has 
generated  a  region  of  turbulence  (Fig.  3.7d),  and  in  the  upper 
part  of  the  old  convective  layer  a  decaying  zone  of  turbulent 
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Fig.  3.6c.  As  Fig.  3.6a,  but  for  the  w-field.  Contour 
interval  is  1  cms-1 . 


i 


illustrated  in  Fig.  3.9  showing  the  fields  of  e  and  i  at  12  h 
in  the  case  where  the  large-scale  wind  is  directed  seawards 
( ug=-3ms_1 ) .  Here  the  turbulence  energy  produced  over  land  de¬ 
cays  as  it  moves  over  the  colder  sea.  The  length  scale  survives 


much  further,  as  shown  in  Fig.  3.9b.  The  mean  flow  in  the  model 
is  influenced  only  by  the  turbulence  fields  e  and  4  through  the 
diffusivity  K  ~*/e,  and  therefore  this  behaviour  of  *  has  only 
a  slight  influence  on  the  mean  flow.  The  importance  of  predicting 
the  turbulence  length  scale  and  the  turbulence  energy  is  primar¬ 
ily  connected  to  the  application  of  the  model  to  air  pollution 
problems.  The  study  of  coastal  meteorology  has  been  concentrated 
on  the  dispersion  conditions  in  many  studies  in  recent  years,  in 
which  case  knowledge  of  both  the  mean  flow  and  the  turbulence 
spectra  are  important.  The  simplest  model  of  turbulence  spectra 
which  can  be  of  use  in  this  case  must  enable  a  determination  of 
a  timescale  (Lagrangian  integral  scale),  and  a  total  variance. 
The  case  of  dispersion  from  a  single  point  in  decaying  turbu- 


90 


Fig.  3.7a.  Contour  plot  of  the  u-component  at  t=  14  h 
(2  hours  after  model  sunset).  Numbers  on  contours  in  ms 
(Reference  case). 


Fig.  3.7b.  As  fig.  3.7a,  but  for  the  v-component. 


the  turbulence 


lence  using  a  simple  two-parameter  model  is  treated  in  Troen  et 
al.  (1982),  where  the  the  relevance  of  the  knowledge  of  both  e 
and  2  is  treated  in  more  detail.  With  essentially  constant 
length  scale  and  decaying  energy  the  timescale  of  the  turbulence 
will  increase,  and  the  model  thus  predicts  a  stage  in  the  decay 
with  weak  and  slow  fluctuations  in  the  turbulence  equivalent  to 
the  meandering  often  observed  under  light  wind  conditions  (Kri- 
stensen  et.  al.  1981).  The  modelling  of  the  length  scale  in  the 
present  model  is  obviously  besed  on  a  rather  arbitrary  prognostic 
equation  which  however,  is  simple  and  has  the  correct  asymptotic 
behaviour  in  transitional  periods,  provided  of  course  that  the 
equilibrium  values  2S  are  correctly  predicted.  Using  instead 
only  the  energy  equation  and  a  diagnostic  value  for  2  equivalent 
to  2=2S  leads  to  a  larger  dissipation  in  the  energy  equation  in 
the  decay  case  considered.  Returning  to  the  behaviour  of  the 
sea-breeze  front  we  present  Fig.  3.8,  which  shows  the  u-field 
at  18  h,  where  the  front  has  reached  approximately  91  km  inland. 
The  maximum  value  for  u  has  further  decreased,  but  the  maximum 
vertical  velocity  (not  shown)  has  decreased  from  18  cms-1  (at 
14  h)  to  only  14  cms-1.  At  this  point  only  very  little  further 
inland  progression  of  the  front  takes  place. 


4.  CONCLUSION 

The  flow  generated  by  periodic  heating  of  the  land  surface  while 
maintaining  a  constant  sea-surface  temperature  (sea-breeze)  has 
been  investigated  using  analytical  and  numerical  modelling.  The 
linearized  equations  of  motion  can  be  nondimensionalized  in  such 
a  way  as  to  make  solutions  universal  except  for  a  dependence  on 
the  scaled  coriolis  parameter  f8»f/»,  where  u>  is  the  frequency 
of  the  surface  forcing.  In  the  case  when  the  breeze  is  superim¬ 
posed  on  a  large  scale  flow  across  the  coast  a  strong  attenuation 
of  the  breeze  circulation  results.  It  is  suggested  that  when  the 
linearization  becomes  invalid,  at  finite  values  of  the  surface 
forcing,  the  effect  on  the  strength  of  the  circulation  can  be 
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taken  into  account  by  a  simple  "self-advection"  model.  The  ana¬ 
lytical  model  has  been  generalized  to  include  the  effects  of  ver¬ 
tical  structure  on  the  diffusivity  and  stratification.  The  effect 
of  an  elevated  inversion  with  a  simultaneous  jump  in  effective 
vertical  mixing  and  stratification  prescribed,  is  a  decrease  of  ! 

the  strength  of  the  circulation  when  the  inversion  height  is  com¬ 
parable  with  the  scale  height.  This  appears  as  a  combined  effect 
of  the  limitation  of  the  vertical  mixing  and  the  increased  strati- 
tification,  which  can  be  studied  independently  in  the  analytical 
model.  The  analytical  method  requires  a  simplification  of  the 
physical  model,  which  in  the  case  of  the  sea-breeze  is  substan¬ 
tial.  The  more  detailed  modelling  of  the  flow  requires  the  use 
of  numerical  techniques,  with  which  the  interplay  of  the  mean 
flow  and  the  turbulent  transport  can  be  handled  more  realisti¬ 
cally.  The  numerical  described  in  the  present  report  employs  a 
turbulence  parameterization  scheme,  based  on  the  turbulence  , 

energy  equation  and  a  simple  prognostic  equation  for  the  turbu¬ 
lence  length  scale.  The  model  is  found  to  be  able  to  reproduce  , 

a  number  of  observed  features  in  sea-breeze  situations,  in  par¬ 
ticular  a  realistic  formation  of  a  sharp  sea-breeze  front  and 
the  inland  propagation  of  this  front.  The  qualitative  response 
of  the  sea-breeze  to  different  synoptic  situations  is  also  found 
to  be  in  qualitative  agreement  with  observations.  The  prediction 
of  the  development  of  both  the  turbulence  energy  and  length  scale 
has  applications  to  dispersion  modelling  in  particular  when  a 
more  accurate  description  than  offered  by  K-models  is  required, 
but  presently  very  little  is  known  about  the  detailed  evolution 
of  the  turbulence  field  under  instationary  conditions,  and  the  >' 

study  presented  here  offers  only  a  simple  model  which  should  be 
tested  by  comparison  with  observational  data  before  definite  con¬ 
clusions  can  be  made  as  to  this  aspect. 

The  increase  in  frontal  speed  relative  to  the  maximum  speed  in 

the  evening  leads  to  a  structure  of  the  flow,  with  a  small  region 

behind  the  progressing  front,  moving  at  a  speed  large  enough  to 

follow  the  front  inland,  without  any  further  air  from  the  sea 

entering  into  this  region.  This  structure  closely  resembles  the 

"cut-off  vortex"  described  by  Simpson  et.  al.  (1977).  A  more  ac-  i 
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curate  modelling  of  this  feature  will  probably  require  a  larger 
resolution  in  the  horizontal. 

An  increase  in  horizontal  resolution  down  to  of  the  order  of 
hundreds  of  meters  is  possible  within  the  nonhydrostatic  model- 
formulation,  with  the  exception  that  the  treatment  of  the  hori¬ 
zontal  diffusion  must  be  improved  using  possibly  a  formulation 
along  the  lines  followed  in  the  treatment  of  the  vertical  dif¬ 
fusion  terms. 


ACKNOWLEDGEMENTS 


The  interest  in  the  analytical  modelling  technique  was  encouraged 
during  a  visit  to  Oregon  State  University,  in  particular  as  a  re¬ 
sult  of  the  support  given  by  Prof.  Larry  Mahrt.  The  visit  was 
made  possible  partially  by  a  grant  given  by  the  Nato  Air-Sea 
Interaction  Panel. 

Drs.  E.W.  Peterson  and  T.  Mikkelsen  made  a  number  of  improvements 
of  the  text  and  helped  in  the  preparation  of  the  manuscript. 


V 


AD  -  A  1 36  702  ANALYTICAL  AND  NUMERICAL  MODELLING  OF  FLOW  ORIVFN  BY 
SURFACE  DIFFERENTIAL  HEAT INGIU)  R I SOE  NATIONAL  l AB 
R0SK1L0E  (OENMARKl  I  TROEN  SEP  82  RISOER  452 

UNCI ASSI F 1  ED  F/G  4/1 


-  97  - 


LIST  OF  SYMBOLS 

a  (fs  +  1 )“1/2. 

Aik  Matrix  defined  by  Eq.  (3.22). 

A(  n)  /A^^  (  n)  Amplitude  function  of  a  u  com¬ 
ponent  corresponding  to  wavenumber  k  and 
frequency  w. 

b  <fs  -  1)"1/2 

Matrix  defined  by  Eq.  (3.22). 

c  Constant  defined  on  p.  19. 


C1  *  c2 »  c3  *  c4 


Constants  in  turbulence  energy  equation 
(p.  66). 


Depth  scale 


Nondimensional  value  of  outer  radius  in 
Fourier-Bessel  expansion  (D  >>  1). 


Divergence  operator  defined  on  p.  73, 


D( u* ,w* ) 


Fourier  component  of  D*  with  wavenumber  k. 


Turbulence  kinetic  energy. 


Equilibrium  turbulence  kinetic  energy  (p.65), 


Coriolis  parameter. 
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Froude  number. 


function  of  (  ) . 


gravitational  acceleration. 


function  of  (  ) . 


Geostrophic  windspeed. 


depth  scale 


H,  Hu 


depth  scale,  *  /K/m 


Height  of  boundary  layer. 


surface  heatflux. 


Height  of  boundary  layer  at  "infinity"  =  xmax . 


imaginary  unit,  index 


Jo*  J 1 


Bessel  functions  of  zero  and  first  order. 


wavenumber  in  x-direction. 


unit  vector  in  z-direction. 


N/U  •  /k7JT. 


Km»  ^h 


dif fusivity 


diffusivity  for  momentum  and  heat. 


mixing  length. 


equilibrium  mixing  length  (p.  65) 


Length  scale  ■  N/»  •  /K/u. 


r  » 

i 
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Monin-Obukhov  length. 

stratification  parameter  (p.  11). 

Brunt-Vaisala  frequency. 

pressure . 

defined  on  p.  10. 

standard  pressure. 

heatf lux  =  <w' 6' > 

heat  source  term. 

scaled  radius. 

constant  =  »^2/2  ( 1  +  i ) 

complex  conjugate  of  r  =  /2/2(1-i). 

Relative  deviation  of  potential  temperature 

»  (9-e)/e 

gas  constant  for  dry  air. 

Bulk  Richardson  number  (p.  67). 
time. 

timescale,  temperature, 
surface  temperature. 


relaxation  timescale  of  surface  temperature 
<p.  74). 

length  of  day  (p.  74). 
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u 


u 


u(C,n,t) 


u(k,n) 

„(*> 


u* 


uj 

U(U) 

umax 


$ 

v(C,n,T) 

v(k,n) 

vj 

v9 


velocity  component  in  x-direction 

dimensional  value  of  u  (in  ch.  2). 

scaled  value  of  u  as  function  of  scaled 
coordinates. 

Fourier  transform  defined  on  p.  9. 

=  Re  {u(k,n)cos  k  Z/k  exp(-iT)}. 
friction  velocity. 

model  u-component  before  adding  3xir-term. 
x- component  of  geostrophic  windspeed. 
defined  on  p.  12-14. 

max  value  of  u  in  the  presence  of  a  large 
advective  velocity  U. 

large-scale  advective  velocity  in  x-direction, 
horizontal  scale  velocity. 

velocity  component  in  y-direction. 

dimensional  value  of  v  (Ch.  2). 

/ 

windvector.  , 

scaled  v  as  function  of  scaled  coordinates. 
Fourier  transform  of  v  defined  on  p.  9. 
defined  on  p.  12-14. 
y- component  of  geostrophic  windvector. 


w 


vertical  component  of  windvector 


w 

w 

w(C,n,x) 
w(k, n) 

w* 

w* 

Wj 

X 


y 

z 

zo 

zi 

a,  aj 

B,Bj 

r 


*c 


rc 

A 


dimensional  value  of  w  (Ch.  2). 
vertical  velocity  scale. 

scaled  w  as  function  of  scaled  coordinates 
Fourier  transform  of  w  defined  on  p.  9. 
convective  velocity  scale  =  (g/T0<w' 6'>H) 1 
model  w- component  before  adding  3 2 *- term, 
defined  on  p.  12-14. 
distance  inland  from  coast, 
distance  along  coast, 
height. 

roughness  length, 
height  of  i'th  gridpoint. 

*  ±  /F  with  Re (a)  <  0. 
characteristic  roots  (p.  10). 

*  (<*j  +  i) 

-  38/3z 

defined  on  p.  66. 
defined  on  p.  66. 
difference  operator. 


/ 
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J 

j 

i 

I 

.  ! 


! 

I 


Az; 


n 

n' 

e 

3 


8U,n,t) 

0(k,n) 

0  j 
e 

x 


grid  spacing  ■  z^+i  -  z^. 

scaled  vertical  coordinate  *  z//k/<d 

scaled  vertical  coordinate  =  z//k/f 

potential  temperature. 

dimensional  value  of  S  (Ch.  2). 

scaled  8  as  function  of  scaled  coordinates. 

Fourier  transform  of  8  defined  on  p.  9. 

defined  on  p.  12-14. 

standard  potential  temperature. 

-  k’/3  n. 


V 


II 


*h 

’k 

e 

4 


♦s 


♦i 

♦ik,,4ik)<n),4<n) 


kinematic  viscosity. 

pressure  function  defined  on  p.  64. 

hydrostatic  pressure  function  (p.  64). 

Fourier  mode  of  wavenumber  k  of  «. 

scaled  x-coordinate. 

latitude. 

defined  on  p.  20. 

"chapeau”  function  defined  on  p.  72. 

phase  function  of  u- component  corresponding 
to  wavenumber  k  and  frequency  «. 


direction  of  wind  maximum  (p.  26). 
density. 

standard  density  (p.  63). 
i'th  root  of  JQ. 
scaled  time  *  wt. 

scaled  time  of  wind  maximum  (p.  26). 
diurnal  frequency,  forcing  frequency 
scaled  diurnal  frequency  *  w/f. 

»  k  U/L  "induced  forcing  frequency" 

diurnal  frequency. 

length  scale  parameter  *  (1/L)2^3. 
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APPENDIX 
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A1 »  Ekman  problem  in  two-layer  model 


In  a  homogeneous  stationary  boundary  layer  the  usual  Ekman  sol¬ 
ution  for  the  horizontal  wind  vector  results  from  the  equations 
of  motion  if  a  height-independent  diffusivity  for  momentum  is 
assumed.  The  solution  for  the  case  where  the  vertical  mixing  is 
limited  to  some  depth  H,  for  example  as  a  consequence  of  an  in¬ 
version  layer  at  this  height,  can  be  obtained  by  assuming  the 
atmosphere  to  be  divided  into  two  layers  with  different  dif- 
fusivities  in  the  two  layers  and  matching  the  solutions  for  each 
layer  at  the  interface.  Using  the  Ekman  depth  corresponding  to 
the  diffusivity  in  the  lowest  layer  for  height  scaling  i.e.  de¬ 
fining  n  =»  z(2K^/f)1/2,  scaling  wind  velocity  by  the  geostrophic 
windspeed  G,  and  using  complex  notation  as  V  »  G-1(u-ug  + 
i(v-vg) )  the  dynamical  equations  car.  be  written: 


v  +  a* 


3ni 


<A1  ) 


where 


1  in  the  lower  layer 

*2 '  { 

K2/K1  in  the  top  layer  . 

The  solution  in  each  layer  can  be  written  as  a  sum  of  two  complex 
exponentials,  and  the  four  coefficients  are  fixed  by  application 
of  the  boundary  conditions  V  *  -1  at  n  *  0  (corresponding  to  u 
and  v  zero  at  the  surface  and  the  alignment  of  the  real  axis 
along  the  direction  of  the  geostrophic  wind  vector),  V  ♦  0  for 
n  ♦  «,  and  the  matching  conditions  at  n  *  b  =  H(2K1/f )_1 ^2  from 
the  requirement  of  continuity  of  the  wind  vector  and  the  momentum 
flux,  viz 


v^b)  =  v2(b) 


(A2 ) 

3Vl  Sv2 

-  =  - —  at  n  =  b 

3  n  3  n 

The  roots  in  the  characteristic  equation  for  Eq.  (A1 )  are 
±  a-1  /TT  and  the  solution  can  be  written 

v  =  A  expU1+i)r>)  =  B  exp(-(1  +  i)n)  (A3) 

for  the  lower  layer,  and  as 

v  =  C  exp(-( 1+i )  1)  ( A  4 ) 

3 

for  the  top  layer. 

The  boundary  and  matching  conditions  results  in  the  following 
equations  for  the  coefficients: 

A  +  B  *  -  1 

Apb  +  Bp~b  =  Cp"b/a  ( A5 ) 

Apb  -  Bp~b  =  -  Cap-b/a 

where  p  *  exp(1+i).  The  solution  to  this  system  of  equations  is 
easily  obtained: 

B  ' '  ft  p2b  (ft  -  ’]  '  <»‘> 

c .  (,  -  ?±\  „(^b)r!ii  . .  r’ 


{ 

!  * 

1 


no 


i 


2 
.  1 

i 

i 

i 
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For  a  *  1  corresponding  to  K2  *  Ki  (and  no  discontinuity)  we 
obtain  A  =  0,  B  =  -1  and  C  *  -1  and  the  usual  Ekman  solution 
given  by  Eq.  (A3).  For  the  case  of  vanishing  vertical  mixing 
in  the  top  layer  we  have  a  *  0,  and  Eq.  (A6)  simplifies  to: 

A  =  -  [p2b  +  1 ]-1 

( A7  ) 

B  =  -  p2b  [p2b  +  1 ]_1 


The  coefficient  C  will  tend  to  infinity  as  a  +0,  but  as  expected 
the  wind  tends  to  the  geostrophic  wind  everywhere  in  the  top 
layer  as  a  ♦  0  as  can  be  seen  by  combining  (A6)  and  (A4).  with 
K2  *  0  undamped  inertial  oscillations  are  also  possible  solutions 
in  the  top  layer.  The  usual  Ekman  solution  is  again  obtained  from 
(A7 )  if  the  discontinuity  tends  to  infinity.  From  (A7)  and  (A3) 
we  finally  obtain  the  wind  profiles  in  the  lower  layer  as 

u/G  *  1  -  D [ cosh ( n-b) cos ( n-b) cosh  b  cosb  + 


sinh  b  sinb  sinh( n-b) sin( n-b) 
v/G  *  D[cosh  v  cosb  sinh( n-b) sin( n-b)  - 
cosh ( n-b) cos ( n-b) sinh  b  sinb] 


(A8) 


with 


D  *  [cosh2b  cos2b  +  sinh2b  sin^bj”1  . 

These  profiles  are  shown  on  Fig.  A1  for  different  values  of  the 
thickness  of  the  mixed  layer  b.  Also  shown  is  the  wind  vector 
at  the  top  of  the  mixed  layer  obtained  from  Eq.  (A8)  by  setting 
n  *  b. 


i 

j 


Fig.  A.1  The  Ekman  spiral  for  different  values  of  the  scaled 
thickness  of  the  boundary  layer  (see  text).  The  numbers  on 
the  curves  give  the  thickness.  The  endpoint  of  each  curve 
corresponds  to  the  top  of  the  boundary  layer,  and  the  dashed 
curve  connects  all  possible  endpoints. 

Geostrophic  value  in  this  (u,v)-plot  is  G»(1,0). 
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